---
title: "Bats and human blood project: gender-specific season analysis"
output:
  html_document:
    df_print: paged
---

```{r setup, echo=FALSE}
knitr::opts_knit$set(root.dir = here::here())
# setwd(here::here())
```

# prepare & load packages
```{r}
# pathWD <- "/Volumes/5360/Biophysics/Experiments/dRT-DC/220301_Bob_DataAnalysis_BatsHumans/20230829_NycNoc2021SexMonthAnalysis"
pathWD <- "~/Documents/08_Pojekte/01_dRT-DC/220301_Bob_DataAnalysis_BatsHumans/20230829_NycNoc2021SexMonthAnalysis"
setwd(pathWD)
# setwd(here::here())

pacman::p_load(tidyverse, lmerTest, car, emmeans, easystats, readr, ggbeeswarm)

# load data ####
Bat1Table <- read_delim(file = "../Bat1Table_matProp.csv",
                        delim = ",",
                        col_names = TRUE)
Bat1Table$replicate <- paste0("NN", Bat1Table$index)

Bat2Table <- read_delim(file = "../Bat2Table_matProp.csv",
                        delim = ",",
                        col_names = TRUE)
Bat2Table$replicate <- paste0("NN", Bat2Table$index)

# concatenate data tables
combTable <- rbind(Bat1Table,Bat2Table)
remove(Bat1Table,Bat2Table)

# only measurement from 2021
combTable <- subset(combTable, combTable$year == 2021)

# combTable
View(combTable)
```

# seasonal effect: month categogries
```{r}
summary(combTable$month)

combTable <- combTable %>% mutate(monthGroup = case_when(month <= 5 ~ "4-5",
                                                         month <= 8 ~ "8",
                                                         month <= 9 ~ "9",
                                                         TRUE ~ "10"))

```

# define factors & sample size tables
```{r}
# define factors ####
combTable$animalType <- as.factor(combTable$animalType)
combTable$replicate <- as.factor(combTable$replicate)
combTable$tempLevel <- as.factor(combTable$tempLevel)
combTable$sex <- as.factor(combTable$sex)
combTable$operator <- as.factor(combTable$operator)
combTable$monthGroup <- factor(combTable$monthGroup,
                               levels = c("4-5","8","9","10")) # sort

# re-level temperature levels: "Warm" "RT" "Cold"
combTable$tempLevel <- relevel(combTable$tempLevel, ref = "RT")
combTable$tempLevel <- relevel(combTable$tempLevel, ref = "Warm")

# levels ####
levels(combTable$animalType)
levels(combTable$replicate)
levels(combTable$tempLevel)
levels(combTable$sex)
levels(combTable$operator)
levels(combTable$monthGroup)

# overview table ####
table(combTable$monthGroup, combTable$replicate)
table(combTable$monthGroup)

table(combTable$monthGroup, combTable$sex)

table(combTable$sex, combTable$replicate)
table(combTable$sex)

```

# model fit

# preparation
```{r}
levels.Month <- levels(combTable$monthGroup)
levels.Temp <- levels(combTable$tempLevel)
levels.Sex <- levels(combTable$sex)

listOutcome <- c("area","tau1","tau2","deltaDefoTau1","deltaDefoTau2","eModulus","viscosity")
```

# density plots
```{r}
plots.Density <- purrr::map(listOutcome,
                         ~ ggplot(data = combTable,
                                  mapping = aes_(x = as.name(.x))) +
                           geom_density(size=1, alpha=0.3) +
                           ggtitle(paste(.x)))
plots.Density
```


# all models
```{r}
# formulas
listFormula.full <- purrr::map(listOutcome,
                               ~ as.formula(paste(.x, "~ monthGroup * tempLevel * sex + (1|replicate)"))
)

# calculate models
listModel.full <- purrr::map(listFormula.full,
                             ~ lmer(formula = .x,
                                    data = combTable,
                                    REML = TRUE,
                                    na.action = na.exclude)
)

# step function from lmerTest package
# remove interactions which are not significant
# but keep always the main effects
listModel.step <- purrr::map(listModel.full,
                             ~ lmerTest::step(.x,
                                              ddf = "Satterthwaite",
                                              alpha.fixed = 0.05,
                                              alpha.random = 0.1,
                                              reduce.fixed = TRUE,
                                              reduce.random = FALSE,
                                              keep = c("monthGroup", "tempLevel", "sex")
                             )
)

# extract final models
listModel.final <- purrr::map(listModel.step,
                              ~ get_model(.x)
)

```


# print outome
```{r}
for (m in 1:length(listModel.full)) {
  writeLines("------------------------------------------------------------")
  writeLines(paste("Outcome variable:", listOutcome[m]))
  # writeLines("------------------------------------------------------------")
  # writeLines(paste("Final model:", listModel.final[[m]]@call, "\n", "\n"))
  
  writeLines("-----Anova3-------------------------------------------------")
  # global p-values
  Anova(listModel.final[[m]], type="3") %>% print()
  
  writeLines("-----Summary------------------------------------------------")
  # summary (estimate) output
  summary(listModel.final[[m]], corr=FALSE) %>% print()
  
  writeLines("-----End----------------------------------------------------\n")
  
}
```

# fitted values
```{r}
listResults.fittedValues <- purrr::map(listModel.final,
                                       ~ data.frame(replicate = combTable$replicate,
                                                    monthGroup = combTable$monthGroup,
                                                    animalType = combTable$animalType,
                                                    sex = combTable$sex,
                                                    tempLevel = combTable$tempLevel,
                                                    fitted = fitted.values(.x)))

# remove NA
listResults.fittedValues <- purrr::map(listResults.fittedValues,
                                       ~ na.omit(.x))

# remove all duplicates
listResults.fittedValues <- purrr::map(listResults.fittedValues,
                                       ~ .x[!duplicated(.x[c('replicate','animalType','tempLevel')]), ])
```


## mean over tempLevel
```{r}
library("dplyr")

listResults.fittedValues.mean <- purrr::map(listResults.fittedValues,
                                            ~ .x %>% group_by(replicate,monthGroup,animalType,sex) %>% summarise_at(vars(fitted), list(fitted.mean = mean)))
```


# plot preparation
```{r}
library("ggh4x")
library("tools") # for capitalization

# capital names for Gender
labels.Sex <- toTitleCase(levels.Sex)
names(labels.Sex) <- levels.Sex

# labels for months
labels.Month <- c("April-May", "August", "September", "October")
names(labels.Month) <- levels.Month
labels.Month.Short <- c("Apr-May", "Aug", "Sept", "Oct")
names(labels.Month.Short) <- levels.Month

# labels for temp
labels.Temp <- c("37", "23", "10")
names(labels.Temp) <- levels.Temp
labels.TempUnit <- paste(labels.Temp, "°C")
names(labels.TempUnit) <- levels.Temp

# axes titles
title.Month <- "Month"
title.Temp <- "Temperature"
title.TempUnit <- paste(title.Temp, "(°C)")
title.Sex <- "Sex"
title.Temp.Vicinity <- "Ambient temp (°C)"

# labels for outcome variables
listLabelOutcome <- c(bquote("Cell size" ~ ("μm"^2)),
                      expression(tau["inlet"] ~ "(ms)"),
                      expression(tau["channel"] ~ "(ms)"),
                      expression(hat(d)["inlet"]),
                      expression(hat(d)["channel"]),
                      "Young's modulus (kPa)",
                      "Viscosity (Pa s)")

```

# month, results
```{r}
listResults.Month <- purrr::map(listModel.final,
                                ~ emmeans(.x,
                                          pairwise ~ monthGroup,
                                          lmer.df = "satterthwaite",
                                          adjust = NULL,
                                          lmerTest.limit = 1e6)
)

listResults.Month.means <- purrr::map(listResults.Month,
                                      ~ .x$emmeans %>% as.data.frame()
)
listResults.Month.means <- purrr::map(listResults.Month.means,
                                       ~ within(.x, {
                                         lower.SEM <- emmean - SE
                                         upper.SEM <- emmean + SE})
)

listResults.Month.contrasts <- purrr::map(listResults.Month,
                                          ~ .x$contrasts %>% as.data.frame()
)

# adjustment for multiple comparisons: Tukey test
# listResults.Month.contrasts <- purrr::map(listResults.Month.contrasts,
#                                           ~ within(.x,
#                                                    p.value.adj <- 1 - ptukey(.x$t.ratio * sqrt(2),
#                                                                              nmeans = nrow(.x),
#                                                                              df = .x$df))
# )
# adjustment for multiple comparisons: Benjamini Hochberg
listResults.Month.contrasts <- purrr::map(listResults.Month.contrasts,
                                          ~ within(.x,
                                                   p.value.adj <- p.adjust(.x$p.value,
                                                                           method = "BH"))
)

rm(listResults.Month)
```


### plot preparation
```{r}
# adapt contrasts data set for plotting
for (m in 1:length(listModel.final)) {
  listResults.Month.contrasts[[m]] <- listResults.Month.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value.adj < 0.001 ~ "***",
                               p.value.adj < 0.01  ~ "**",
                               p.value.adj < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value.adj)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = c(1,1,1,2,2,3) + 0.05,
           xEnd     = c(2,3,4,3,4,4) - 0.05,
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots month v2
```{r fig.height=8/2.54, fig.width=9.3/2.54}
plots.Month.limits <- data.frame(min = c(17, 2, 2, 0.07, 0.07, 0.90, 1.4),
                                 max = c(24.6, 11, 11, 0.25, 0.25, 1.10, 3.2))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.Month.contrasts[[m]] <- listResults.Month.contrasts[[m]] %>%
    mutate(yStart = (0.69 + 0.08*c(0,1,3,0,2,0)) * (plots.Month.limits$max[m] - plots.Month.limits$min[m]) + plots.Month.limits$min[m],
           yEnd = yStart,
           yText = yStart + 0.015*(plots.Month.limits$max[m] - plots.Month.limits$min[m]))
}

plots.Month <- list()
tbl.plotResults.Month <- data.frame(matrix(ncol = 4, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.Month[[m]] <- ggplot(data = listResults.Month.means[[m]])
  
  # define layout parameters
  plots.Month[[m]] <- plots.Month[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")
  
  # create point plot: replicate means
  plots.Month[[m]] <- plots.Month[[m]] +
    geom_beeswarm(data = listResults.fittedValues.mean[[m]],
                  mapping = aes(x = monthGroup,
                                y = fitted.mean),
                  colour = "#C0C0C0",
                  size = 2,
                  cex = 1.6)
  # are some means out of viewing window?
  if( (plots.Month.limits$min[m] > min(listResults.fittedValues.mean[[m]]$fitted.mean)) | (plots.Month.limits$max[m] < max(listResults.fittedValues.mean[[m]]$fitted.mean)) ) {
    print(paste("Replicate means out of view window:",listLabelOutcome[m]))
  }
  
  # create LMM plot: mean
  plots.Month[[m]] <- plots.Month[[m]] +
    geom_point(mapping = aes(x = monthGroup,
                             y = emmean),
               stat = "identity", # important since we provide mean values already
               size = 3,
               colour = "#1f77bd")
  
  # add LMM plot: error bars
  plots.Month[[m]] <- plots.Month[[m]] +
    geom_errorbar(mapping = aes(x = monthGroup,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = "#1f77bd",
                  position = position_dodge(0.1))
  
  # add discrete x-axis labels
  plots.Month[[m]] <- plots.Month[[m]] +
    scale_x_discrete(labels = labels.Month.Short)
  
  # add p-values
  plots.Month[[m]] <- plots.Month[[m]] +
    geom_segment(data = listResults.Month.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.Month.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.Month[[m]] <- plots.Month[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.Month[[m]] <- plots.Month[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.Month.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Month,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.Month[[m]] <- plots.Month[[m]] +
    coord_cartesian(ylim = c(plots.Month.limits$min[m], plots.Month.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Month", ".svg")),
         width = 9.3, height = 8, units = "cm", dpi = 600, device = "svg",
         plot = plots.Month[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Month", ".png")),
         width = 9.3, height = 8, units = "cm", dpi = 600, device = "png",
         plot = plots.Month[[m]], scale = 1)
  
  # results table
  tbl.plotResults.Month <- rbind(tbl.plotResults.Month, data.frame(rep(listOutcome[m], times = nrow(listResults.Month.means[[m]])), listResults.Month.means[[m]]$monthGroup, listResults.Month.means[[m]]$emmean, listResults.Month.means[[m]]$SE, listResults.Month.means[[m]]$upper.CL - listResults.Month.means[[m]]$emmean))
}
 
colnames(tbl.plotResults.Month) <- c("targetVar", "monthGroup", "mean", "SEM", "CI")
write.csv(tbl.plotResults.Month, file.path(pathWD, paste0("tblPlotResults", "_Month", ".csv")), row.names = FALSE)

plots.Month

```


# sex, results
```{r}
listResults.Sex <- purrr::map(listModel.final,
                              ~ emmeans(.x,
                                        pairwise ~ sex,
                                        lmer.df = "satterthwaite",
                                        adjust = NULL,
                                        lmerTest.limit = 1e6)
)

listResults.Sex.means <- purrr::map(listResults.Sex,
                                    ~ .x$emmeans %>% as.data.frame()
)
listResults.Sex.means <- purrr::map(listResults.Sex.means,
                                     ~ within(.x, {
                                       lower.SEM <- emmean - SE
                                       upper.SEM <- emmean + SE})
)

listResults.Sex.contrasts <- purrr::map(listResults.Sex,
                                        ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.Sex)
```

### plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.Sex.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listResults.Sex.contrasts[[m]] <- listResults.Sex.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = c(1) + 0.05,
           xEnd     = c(2) - 0.05,
           xText  = 0.5*(xStart+xEnd))
}
```

## all plots sex
```{r}
plots.Sex.limits <- data.frame(min = c(15, 2, 2, 0.07, 0.07, 0.85, 1.4),
                                max = c(25, 11, 11, 0.25, 0.25, 1.05, 3.2))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.Sex.contrasts[[m]] <- listResults.Sex.contrasts[[m]] %>%
    mutate(yStart = rep((0.85 + 0.08*c(1)) * (plots.Sex.limits$max[m] - plots.Sex.limits$min[m]) + plots.Sex.limits$min[m], times = 1),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.Sex.limits$max[m] - plots.Sex.limits$min[m]))
}

plots.Sex <- list()
tbl.plotResults.Sex <- data.frame(matrix(ncol = 4, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.Sex[[m]] <- ggplot(data = listResults.Sex.means[[m]])
  
  # define layout parameters
  plots.Sex[[m]] <- plots.Sex[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")
  
  # create plot
  plots.Sex[[m]] <- plots.Sex[[m]] +
    geom_bar(mapping = aes(x = sex,
                           y = emmean,
                           fill = sex),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.Sex[[m]] <- plots.Sex[[m]] +
    geom_errorbar(mapping = aes(x = sex,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.Sex[[m]] <- plots.Sex[[m]] +
    scale_x_discrete(labels = labels.Sex)
  
  # add p-values
  plots.Sex[[m]] <- plots.Sex[[m]] +
    geom_segment(data = listResults.Sex.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.Sex.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.Sex[[m]] <- plots.Sex[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.Sex[[m]] <- plots.Sex[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.Sex.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Sex,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.Sex[[m]] <- plots.Sex[[m]] +
    coord_cartesian(ylim = c(plots.Sex.limits$min[m], plots.Sex.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Sex", ".svg")),
         width = 6, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.Sex[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Sex", ".png")),
         width = 6, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.Sex[[m]], scale = 1)
  
  # results table
  tbl.plotResults.Sex <- rbind(tbl.plotResults.Sex, data.frame(rep(listOutcome[m], times = nrow(listResults.Sex.means[[m]])), listResults.Sex.means[[m]]$sex, listResults.Sex.means[[m]]$emmean, listResults.Sex.means[[m]]$SE, listResults.Sex.means[[m]]$upper.CL - listResults.Sex.means[[m]]$emmean))
}

colnames(tbl.plotResults.Sex) <- c("targetVar", "sex", "mean", "SEM", "CI")
write.csv(tbl.plotResults.Sex, file.path(pathWD, paste0("tblPlotResults", "_Sex", ".csv")), row.names = FALSE)

plots.Sex

```


# temp, results
```{r}
listResults.Temp <- purrr::map(listModel.final,
                               ~ emmeans(.x,
                                         pairwise ~ tempLevel,
                                         lmer.df = "satterthwaite",
                                         adjust = NULL,
                                         lmerTest.limit = 1e6)
)

listResults.Temp.means <- purrr::map(listResults.Temp,
                                     ~ .x$emmeans %>% as.data.frame()
)
listResults.Temp.means <- purrr::map(listResults.Temp.means,
                                     ~ within(.x, {
                                       lower.SEM <- emmean - SE
                                       upper.SEM <- emmean + SE})
)

listResults.Temp.contrasts <- purrr::map(listResults.Temp,
                                         ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.Temp)
```


### plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.Temp.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listResults.Temp.contrasts[[m]] <- listResults.Temp.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = c(1,1,2) + 0.05,
           xEnd     = c(2,3,3) - 0.05,
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots temp
```{r}
plots.Temp.limits <- data.frame(min = c(15, 2, 2, 0.07, 0.07, 0.5, 0.5),
                                max = c(25, 11, 11, 0.25, 0.25, 1.6, 3.8))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.Temp.contrasts[[m]] <- listResults.Temp.contrasts[[m]] %>%
    mutate(yStart = rep((0.85 + 0.08*c(0,1,0)) * (plots.Temp.limits$max[m] - plots.Temp.limits$min[m]) + plots.Temp.limits$min[m], times = 1),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.Temp.limits$max[m] - plots.Temp.limits$min[m]))
}

plots.Temp <- list()
tbl.plotResults.Temp <- data.frame(matrix(ncol = 4, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.Temp[[m]] <- ggplot(data = listResults.Temp.means[[m]])
  
  # define layout parameters
  plots.Temp[[m]] <- plots.Temp[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")
  
  # create plot
  plots.Temp[[m]] <- plots.Temp[[m]] +
    geom_bar(mapping = aes(x = tempLevel,
                           y = emmean,
                           fill = tempLevel),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.Temp[[m]] <- plots.Temp[[m]] +
    geom_errorbar(mapping = aes(x = tempLevel,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.Temp[[m]] <- plots.Temp[[m]] +
    scale_x_discrete(labels = labels.Temp)
  
  # add p-values
  plots.Temp[[m]] <- plots.Temp[[m]] +
    geom_segment(data = listResults.Temp.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.Temp.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.Temp[[m]] <- plots.Temp[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.Temp[[m]] <- plots.Temp[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.Temp.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.TempUnit,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.Temp[[m]] <- plots.Temp[[m]] +
    coord_cartesian(ylim = c(plots.Temp.limits$min[m], plots.Temp.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Temp", ".svg")),
         width = 7.5, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.Temp[[m]], scale = 1)
  
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Temp", ".png")),
         width = 7.5, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.Temp[[m]], scale = 1)
  
  # results table
  tbl.plotResults.Temp <- rbind(tbl.plotResults.Temp, data.frame(rep(listOutcome[m], times = nrow(listResults.Temp.means[[m]])), listResults.Temp.means[[m]]$tempLevel, listResults.Temp.means[[m]]$emmean, listResults.Temp.means[[m]]$SE, listResults.Temp.means[[m]]$upper.CL - listResults.Temp.means[[m]]$emmean))
}

colnames(tbl.plotResults.Temp) <- c("targetVar", "tempLevel", "mean", "SEM", "CI")
write.csv(tbl.plotResults.Temp, file.path(pathWD, paste0("tblPlotResults", "_Temp", ".csv")), row.names = FALSE)

plots.Temp

```


# month | sex, results
```{r}
listResults.MonthSex <- purrr::map(listModel.final,
                                   ~ emmeans(.x,
                                             pairwise ~ monthGroup | sex,
                                             lmer.df = "satterthwaite",
                                             adjust = NULL,
                                             lmerTest.limit = 1e6)
)

listResults.MonthSex.means <- purrr::map(listResults.MonthSex,
                                         ~ .x$emmeans %>% as.data.frame()
)
listResults.MonthSex.means <- purrr::map(listResults.MonthSex.means,
                                         ~ within(.x, {
                                           lower.SEM <- emmean - SE
                                           upper.SEM <- emmean + SE})
)

listResults.MonthSex.contrasts <- purrr::map(listResults.MonthSex,
                                             ~ .x$contrasts %>% as.data.frame()
)

# adjustment for multiple comparisons: Tukey test
listResults.MonthSex.contrasts <- purrr::map(listResults.MonthSex.contrasts,
                                             ~ within(.x,
                                                      p.value.adj <- 1 - ptukey(.x$t.ratio * sqrt(2),
                                                                                nmeans = nrow(.x) / length(levels.Sex),
                                                                                df = .x$df))
)
```

### plot preparation
```{r}
yMinFrac <- c(0.7, 0.1, 0.1, 0.5, 0.5, 0.84, 0.4)
sigFrac <- 0.4

# adapt contrasts data set for plotting
listResults.MonthSex.plotPara <- list()
n = nrow(listResults.MonthSex.contrasts[[1]])
for (m in 1:length(listModel.final)) {

  dataMax <- listResults.MonthSex.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.MonthSex.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo     = yLo,
                                               yHi     = yHi)
  
  listResults.MonthSex.contrasts[[m]] <- listResults.MonthSex.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value.adj < 0.001 ~ "***",
                               p.value.adj < 0.01  ~ "**",
                               p.value.adj < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value.adj)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1,1,1,2,2,3) + 0.05, times = length(levels.Sex)),
           xEnd     = rep(c(2,3,4,3,4,4) - 0.05, times = length(levels.Sex)),
           yLineRel = rep((1 + c(0,1,3,0,2,0)) / 4.8, times = length(levels.Sex)),
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi)
  
}
```


## all plots month | sex
```{r}
plots.MonthSex <- list()
tbl.plotResults.MonthSex <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  dataMax <- listResults.MonthSex.plotPara[[m]]$dataMax
  yLo <- listResults.MonthSex.plotPara[[m]]$yLo
  yHi <- listResults.MonthSex.plotPara[[m]]$yHi
  
  plots.MonthSex[[m]] <- ggplot(data = listResults.MonthSex.means[[m]])
  
  # define layout parameters
  plots.MonthSex[[m]] <- plots.MonthSex[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
    facet_grid(cols = vars(sex), labeller = labeller(sex = labels.Sex))
  
  # create plot
  plots.MonthSex[[m]] <- plots.MonthSex[[m]] +
    geom_bar(mapping = aes(x = monthGroup,
                           y = emmean,
                           fill = monthGroup),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.MonthSex[[m]] <- plots.MonthSex[[m]] +
    geom_errorbar(mapping = aes(x = monthGroup,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  # plots.MonthSex[[m]] <- plots.MonthSex[[m]] +
  #   scale_x_discrete(guide = "axis_nested") # guide type comes from package "ggh4x"
  
  # add p-values
  plots.MonthSex[[m]] <- plots.MonthSex[[m]] +
    geom_segment(data = listResults.MonthSex.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.MonthSex.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.MonthSex[[m]] <- plots.MonthSex[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.MonthSex[[m]] <- plots.MonthSex[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.Sex,
         caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Month,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.MonthSex[[m]] <- plots.MonthSex[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_MonthSex", ".svg")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.MonthSex[[m]], scale = 1)
  
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_MonthSex", ".png")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.MonthSex[[m]], scale = 1)
  
  # results table
  tbl.plotResults.MonthSex <- rbind(tbl.plotResults.MonthSex, data.frame(rep(listOutcome[m], times = nrow(listResults.MonthSex.means[[1]])), listResults.MonthSex.means[[m]]$monthGroup, listResults.MonthSex.means[[m]]$sex, listResults.MonthSex.means[[m]]$emmean, listResults.MonthSex.means[[m]]$SE, listResults.MonthSex.means[[m]]$upper.CL - listResults.MonthSex.means[[m]]$emmean))
}

colnames(tbl.plotResults.MonthSex) <- c("targetVar", "monthGroup", "sex", "mean", "SEM", "CI")
write.csv(tbl.plotResults.MonthSex, file.path(pathWD, paste0("tblPlotResults", "_MonthSex", ".csv")), row.names = FALSE)

plots.MonthSex

```


# month | temp, results
```{r}
listResults.MonthTemp <- purrr::map(listModel.final,
                                    ~ emmeans(.x,
                                             pairwise ~ monthGroup | tempLevel,
                                             lmer.df = "satterthwaite",
                                             adjust = NULL,
                                             lmerTest.limit = 1e6)
)

listResults.MonthTemp.means <- purrr::map(listResults.MonthTemp,
                                          ~ .x$emmeans %>% as.data.frame()
)
listResults.MonthTemp.means <- purrr::map(listResults.MonthTemp.means,
                                          ~ within(.x, {
                                            lower.SEM <- emmean - SE
                                            upper.SEM <- emmean + SE})
)

listResults.MonthTemp.contrasts <- purrr::map(listResults.MonthTemp,
                                              ~ .x$contrasts %>% as.data.frame()
)

# adjustment for multiple comparisons: Tukey test
# listResults.MonthTemp.contrasts <- purrr::map(listResults.MonthTemp.contrasts,
#                                              ~ within(.x,
#                                                       p.value.adj <- 1 - ptukey(.x$t.ratio * sqrt(2),
#                                                                                 nmeans = nrow(.x) / length(levels.Temp),
#                                                                                 df = .x$df))
# )
# adjustment for multiple comparisons: Benjamini Hochberg
listResults.MonthTemp.contrasts <- purrr::map(listResults.MonthTemp.contrasts,
                                              ~ within(.x,
                                                       p.value.adj <- p.adjust(.x$p.value,
                                                                               method = "BH"))
)

rm(listResults.MonthTemp)
```


### plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.MonthTemp.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listResults.MonthTemp.contrasts[[m]] <- listResults.MonthTemp.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value.adj < 0.001 ~ "***",
                               p.value.adj < 0.01  ~ "**",
                               p.value.adj < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value.adj)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1,1,1,2,2,3) + 0.05, times = length(levels.Temp)),
           xEnd     = rep(c(2,3,4,3,4,4) - 0.05, times = length(levels.Temp)),
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots month | temp
```{r fig.height=16/2.54, fig.width=10/2.54}
plots.MonthTemp.limits <- data.frame(min = c(15, 1, 1, 0.07, 0.07, 0.5, 0.5),
                                     max = c(26, 12, 12, 0.25, 0.25, 1.85, 5.5))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.MonthTemp.contrasts[[m]] <- listResults.MonthTemp.contrasts[[m]] %>%
    mutate(yStart = rep((0.69 + 0.08*c(0,1,3,0,2,0)) * (plots.MonthTemp.limits$max[m] - plots.MonthTemp.limits$min[m]) + plots.MonthTemp.limits$min[m], times = length(levels.Temp)),
           yEnd = yStart,
           yText = yStart + 0.002*(plots.MonthTemp.limits$max[m] - plots.MonthTemp.limits$min[m]))
}

plots.MonthTemp <- list()
tbl.plotResults.MonthTemp <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.MonthTemp[[m]] <- ggplot(data = listResults.MonthTemp.means[[m]])
  
  # define layout parameters
  plots.MonthTemp[[m]] <- plots.MonthTemp[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 11, angle = 90),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
    facet_grid(rows = vars(tempLevel), labeller = labeller(tempLevel = labels.TempUnit))
  
  # create point plot: replicate means
  plots.MonthTemp[[m]] <- plots.MonthTemp[[m]] +
    geom_beeswarm(data = listResults.fittedValues[[m]],
                  mapping = aes(x = monthGroup,
                                y = fitted),
                  colour = "#C0C0C0",
                  size = 2,
                  cex = 1.6)
  # are some means out of viewing window?
  if( (plots.MonthTemp.limits$min[m] > min(listResults.fittedValues[[m]]$fitted)) | (plots.MonthTemp.limits$max[m] < max(listResults.fittedValues[[m]]$fitted)) ) {
    print(paste("Replicate means out of view window:",listLabelOutcome[m]))
  }
  
  # create LMM plot: mean
  plots.MonthTemp[[m]] <- plots.MonthTemp[[m]] +
    geom_point(mapping = aes(x = monthGroup,
                             y = emmean),
             stat = "identity", # important since we provide mean values already
             size = 3, # line width of bar borders
             colour = "#1f77bd") # bar border color
  
  # add LMM plot: error bars
  plots.MonthTemp[[m]] <- plots.MonthTemp[[m]] +
    geom_errorbar(mapping = aes(x = monthGroup,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = "#1f77bd",
                  position = position_dodge(0.1)) 
  
  # add discrete x-axis labels
  plots.MonthTemp[[m]] <- plots.MonthTemp[[m]] +
    scale_x_discrete(labels = labels.Month.Short)
  
  # add p-values
  plots.MonthTemp[[m]] <- plots.MonthTemp[[m]] +
    geom_segment(data = listResults.MonthTemp.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.MonthTemp.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.MonthTemp[[m]] <- plots.MonthTemp[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.MonthTemp[[m]] <- plots.MonthTemp[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.MonthTemp.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Month,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.MonthTemp[[m]] <- plots.MonthTemp[[m]] +
    coord_cartesian(ylim = c(plots.MonthTemp.limits$min[m], plots.MonthTemp.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_MonthTemp", ".svg")),
         width = 10, height = 16, units = "cm", dpi = 600, device = "svg",
         plot = plots.MonthTemp[[m]], scale = 1)
  
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_MonthTemp", ".png")),
         width = 10, height = 16, units = "cm", dpi = 600, device = "png",
         plot = plots.MonthTemp[[m]], scale = 1)
  
  # results table
  tbl.plotResults.MonthTemp <- rbind(tbl.plotResults.MonthTemp, data.frame(rep(listOutcome[m], times = nrow(listResults.MonthTemp.means[[1]])), listResults.MonthTemp.means[[m]]$monthGroup, listResults.MonthTemp.means[[m]]$tempLevel, listResults.MonthTemp.means[[m]]$emmean, listResults.MonthTemp.means[[m]]$SE, listResults.MonthTemp.means[[m]]$upper.CL - listResults.MonthTemp.means[[m]]$emmean))
}

colnames(tbl.plotResults.MonthTemp) <- c("targetVar", "monthGroup", "tempLevel", "mean", "SEM", "CI")
write.csv(tbl.plotResults.MonthTemp, file.path(pathWD, paste0("tblPlotResults", "_MonthTemp", ".csv")), row.names = FALSE)

plots.MonthTemp

```


# sex | month, results
```{r}
listResults.SexMonth <- purrr::map(listModel.final,
                                   ~ emmeans(.x,
                                             pairwise ~ sex | monthGroup,
                                             lmer.df = "satterthwaite",
                                             adjust = NULL,
                                             lmerTest.limit = 1e6)
)

listResults.SexMonth.means <- purrr::map(listResults.SexMonth,
                                         ~ .x$emmeans %>% as.data.frame()
)
listResults.SexMonth.means <- purrr::map(listResults.SexMonth.means,
                                         ~ within(.x, {
                                           lower.SEM <- emmean - SE
                                           upper.SEM <- emmean + SE})
)

listResults.SexMonth.contrasts <- purrr::map(listResults.SexMonth,
                                             ~ .x$contrasts %>% as.data.frame()
)
```

### plot preparation
```{r}
yMinFrac <- c(0.4, 0.1, 0.1, 0.4, 0.4, 0.3, 0.1)
sigFrac <- 0.2

# adapt contrasts data set for plotting
listResults.SexMonth.plotPara <- list()
n = nrow(listResults.SexMonth.contrasts[[1]])
for (m in 1:length(listModel.final)) {

  dataMax <- listResults.SexMonth.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.SexMonth.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo     = yLo,
                                               yHi     = yHi)
  
  listResults.SexMonth.contrasts[[m]] <- listResults.SexMonth.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1) + 0.05, times = length(levels.Month)),
           xEnd     = rep(c(2) - 0.05, times = length(levels.Month)),
           yLineRel = rep((1 + c(0) ) / 1.8, times = length(levels.Month)),
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi)

}
```


## all plots sex | month
```{r}
plots.SexMonth <- list()
tbl.plotResults.SexMonth <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  dataMax <- listResults.SexMonth.plotPara[[m]]$dataMax
  yLo <- listResults.SexMonth.plotPara[[m]]$yLo
  yHi <- listResults.SexMonth.plotPara[[m]]$yHi
  
  plots.SexMonth[[m]] <- ggplot(data = listResults.SexMonth.means[[m]])
  
  # define layout parameters
  plots.SexMonth[[m]] <- plots.SexMonth[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
    facet_grid(cols = vars(monthGroup))
  
  # create plot
  plots.SexMonth[[m]] <- plots.SexMonth[[m]] +
    geom_bar(mapping = aes(x = sex,
                           y = emmean,
                           fill = sex),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.SexMonth[[m]] <- plots.SexMonth[[m]] +
    geom_errorbar(mapping = aes(x = sex,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.SexMonth[[m]] <- plots.SexMonth[[m]] +
    scale_x_discrete(labels = labels.Sex)
  
  # add p-values
  plots.SexMonth[[m]] <- plots.SexMonth[[m]] +
    geom_segment(data = listResults.SexMonth.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.SexMonth.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.SexMonth[[m]] <- plots.SexMonth[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.SexMonth[[m]] <- plots.SexMonth[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.Month,
         caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Sex,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.SexMonth[[m]] <- plots.SexMonth[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_SexMonth", ".svg")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.SexMonth[[m]], scale = 1)
  
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_SexMonth", ".png")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.SexMonth[[m]], scale = 1)
  
  # results table
  tbl.plotResults.SexMonth <- rbind(tbl.plotResults.SexMonth, data.frame(rep(listOutcome[m], times = nrow(listResults.SexMonth.means[[1]])), listResults.SexMonth.means[[m]]$sex, listResults.SexMonth.means[[m]]$monthGroup, listResults.SexMonth.means[[m]]$emmean, listResults.SexMonth.means[[m]]$SE, listResults.SexMonth.means[[m]]$upper.CL - listResults.SexMonth.means[[m]]$emmean))
}

colnames(tbl.plotResults.SexMonth) <- c("targetVar", "sex", "monthGroup", "mean", "SEM", "CI")
write.csv(tbl.plotResults.SexMonth, file.path(pathWD, paste0("tblPlotResults", "_SexMonth", ".csv")), row.names = FALSE)

plots.SexMonth

```


# temp | month, results
```{r}
listResults.TempMonth <- purrr::map(listModel.final,
                                       ~ emmeans(.x,
                                                 pairwise ~ tempLevel | monthGroup,
                                                 lmer.df = "satterthwaite",
                                                 adjust = NULL,
                                                 lmerTest.limit = 1e6)
)

listResults.TempMonth.means <- purrr::map(listResults.TempMonth,
                                             ~ .x$emmeans %>% as.data.frame()
)
listResults.TempMonth.means <- purrr::map(listResults.TempMonth.means,
                                             ~ within(.x, {
                                               lower.SEM <- emmean - SE
                                               upper.SEM <- emmean + SE})
)

listResults.TempMonth.contrasts <- purrr::map(listResults.TempMonth,
                                                 ~ .x$contrasts %>% as.data.frame()
)
```

### plot preparation (facet grid)
```{r}
yMinFrac <- c(0.8, 0.4, 0.4, 0.6, 0.6, 0.3, 0.2)
sigFrac <- 0.2

# adapt contrasts data set for plotting
listResults.TempMonth_fg.plotPara <- list()
listResults.TempMonth_fg.contrasts <- list()
n = nrow(listResults.TempMonth.contrasts[[1]])
for (m in 1:length(listModel.final)) {

  dataMax <- listResults.TempMonth.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.TempMonth_fg.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo     = yLo,
                                               yHi     = yHi)
  
  listResults.TempMonth_fg.contrasts[[m]] <- listResults.TempMonth.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1,1,2) + 0.05, times = length(levels.Month)),
           xEnd     = rep(c(2,3,3) - 0.05, times = length(levels.Month)),
           yLineRel = rep((1 + c(0,1,0) ) / 2.8, times = length(levels.Month)),
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi)

}
```


## all plots temp | month (facet grid)
```{r}
plots.TempMonth_fg <- list()
tbl.plotResults.TempMonth <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  dataMax <- listResults.TempMonth_fg.plotPara[[m]]$dataMax
  yLo <- listResults.TempMonth_fg.plotPara[[m]]$yLo
  yHi <- listResults.TempMonth_fg.plotPara[[m]]$yHi
  
  plots.TempMonth_fg[[m]] <- ggplot(data = listResults.TempMonth.means[[m]])
  
  # define layout parameters
  plots.TempMonth_fg[[m]] <- plots.TempMonth_fg[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
  facet_grid(cols = vars(monthGroup))
  
  # create plot
  plots.TempMonth_fg[[m]] <- plots.TempMonth_fg[[m]] +
    geom_bar(mapping = aes(x = tempLevel,
                           y = emmean,
                           fill = tempLevel),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.TempMonth_fg[[m]] <- plots.TempMonth_fg[[m]] +
    geom_errorbar(mapping = aes(x = tempLevel,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())

  # add nested x-axis
  # plots.TempMonth_fg[[m]] <- plots.TempMonth_fg[[m]] +
  #   scale_x_discrete(guide = "axis_nested") # guide type comes from package "ggh4x"

  # add p-values
  plots.TempMonth_fg[[m]] <- plots.TempMonth_fg[[m]] +
    geom_segment(data = listResults.TempMonth_fg.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd),
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.TempMonth_fg.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)

  # add nicer fill color
  plots.TempMonth_fg[[m]] <- plots.TempMonth_fg[[m]] +
    scale_fill_brewer(palette = "Blues")

  # add further layout parameter
  plots.TempMonth_fg[[m]] <- plots.TempMonth_fg[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.Month,
       caption = paste("Model:", format(listFormula.full[[m]])),
       x = title.Temp,
       y = listLabelOutcome[m])

  # change coordinates = "zoom in"
  plots.TempMonth_fg[[m]] <- plots.TempMonth_fg[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))

  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempMonth_fg", ".svg")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.TempMonth_fg[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempMonth_fg", ".png")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.TempMonth_fg[[m]], scale = 1)

  # results table
  tbl.plotResults.TempMonth <- rbind(tbl.plotResults.TempMonth, data.frame(rep(listOutcome[m], times = nrow(listResults.TempMonth.means[[1]])), listResults.TempMonth.means[[m]]$tempLevel, listResults.TempMonth.means[[m]]$monthGroup, listResults.TempMonth.means[[m]]$emmean, listResults.TempMonth.means[[m]]$SE, listResults.TempMonth.means[[m]]$upper.CL - listResults.TempMonth.means[[m]]$emmean))
}

colnames(tbl.plotResults.TempMonth) <- c("targetVar", "tempLevel", "monthGroup", "mean", "SEM", "CI")
write.csv(tbl.plotResults.TempMonth, file.path(pathWD, paste0("tblPlotResults", "_TempMonth", ".csv")), row.names = FALSE)

plots.TempMonth_fg

```


### plot preparation (interaction)
```{r}
yMinFrac <- c(0.8, 0.4, 0.4, 0.6, 0.6, 0.3, 0.2)
sigFrac <- 0.2

# adapt contrasts data set for plotting
listResults.TempMonth.plotPara <- list()
n = nrow(listResults.TempMonth.contrasts[[1]])
for (m in 1:length(listModel.final)) {

  dataMax <- listResults.TempMonth.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.TempMonth.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo     = yLo,
                                               yHi     = yHi)
  
  listResults.TempMonth.contrasts[[m]] <- listResults.TempMonth.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = c(1,1,2) + 0.05 + 3*rep(0:(length(levels.Month)-1), each = 3),
           xEnd     = c(2,3,3) - 0.05 + 3*rep(0:(length(levels.Month)-1), each = 3),
           yLineRel = rep((1 + c(0,1,0) ) / 2.8, times = length(levels.Month)),
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi)

}
```


## all plots temp | month (interaction)
```{r}
plots.TempMonth <- list()

for (m in 1:length(listModel.final)) {
  
  dataMax <- listResults.TempMonth.plotPara[[m]]$dataMax
  yLo <- listResults.TempMonth.plotPara[[m]]$yLo
  yHi <- listResults.TempMonth.plotPara[[m]]$yHi
  
  plots.TempMonth[[m]] <- ggplot(data = listResults.TempMonth.means[[m]])
  
  # define layout parameters
  plots.TempMonth[[m]] <- plots.TempMonth[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")
  
  # create plot
  plots.TempMonth[[m]] <- plots.TempMonth[[m]] +
    geom_bar(mapping = aes(x = interaction(tempLevel, monthGroup),
                           y = emmean,
                           fill = tempLevel),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.TempMonth[[m]] <- plots.TempMonth[[m]] +
    geom_errorbar(mapping = aes(x = interaction(tempLevel, monthGroup),
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())

  # add nested x-axis
  plots.TempMonth[[m]] <- plots.TempMonth[[m]] +
    scale_x_discrete(guide = "axis_nested") # guide type comes from package "ggh4x"

  # add p-values
  plots.TempMonth[[m]] <- plots.TempMonth[[m]] +
    geom_segment(data = listResults.TempMonth.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd),
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.TempMonth.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)

  # add nicer fill color
  plots.TempMonth[[m]] <- plots.TempMonth[[m]] +
    scale_fill_brewer(palette = "Blues")

  # add further layout parameter
  plots.TempMonth[[m]] <- plots.TempMonth[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.Month,
       caption = paste("Model:", format(listFormula.full[[m]])),
       x = title.Temp,
       y = listLabelOutcome[m])

  # change coordinates = "zoom in"
  plots.TempMonth[[m]] <- plots.TempMonth[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))

  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempMonth", ".svg")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.TempMonth[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempMonth", ".png")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.TempMonth[[m]], scale = 1)
}

plots.TempMonth

```


# month | temp * sex, results
```{r}
listResults.MonthTempSex <- purrr::map(listModel.final,
                                       ~ emmeans(.x,
                                                 pairwise ~ monthGroup | tempLevel * sex,
                                                 lmer.df = "satterthwaite",
                                                 adjust = NULL,
                                                 lmerTest.limit = 1e6)
)

listResults.MonthTempSex.means <- purrr::map(listResults.MonthTempSex,
                                             ~ .x$emmeans %>% as.data.frame()
)
listResults.MonthTempSex.means <- purrr::map(listResults.MonthTempSex.means,
                                             ~ within(.x, {
                                               lower.SEM <- emmean - SE
                                               upper.SEM <- emmean + SE})
)

listResults.MonthTempSex.contrasts <- purrr::map(listResults.MonthTempSex,
                                                 ~ .x$contrasts %>% as.data.frame()
)

# adjustment for multiple comparisons: Tukey test
listResults.MonthTempSex.contrasts <- purrr::map(listResults.MonthTempSex.contrasts,
                                                 ~ within(.x,
                                                          p.value.adj <- 1 - ptukey(.x$t.ratio * sqrt(2),
                                                                                    nmeans = nrow(.x) / length(levels.Temp) / length(levels.Sex),
                                                                                    df = .x$df))
)
```


### plot preparation
```{r}
yMinFrac <- c(0.6, 0.1, 0.1, 0.4, 0.4, 0.3, 0.2)
sigFrac <- 0.4

# adapt contrasts data set for plotting
listResults.MonthTempSex.plotPara <- list()
n = nrow(listResults.MonthTempSex.contrasts[[1]])
for (m in 1:length(listModel.final)) {

  dataMax <- listResults.MonthTempSex.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.MonthTempSex.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo     = yLo,
                                               yHi     = yHi)
  
  listResults.MonthTempSex.contrasts[[m]] <- listResults.MonthTempSex.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value.adj < 0.001 ~ "***",
                               p.value.adj < 0.01  ~ "**",
                               p.value.adj < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value.adj)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart    = rep(c(1,1,1,2,2,3) + 0.05, times = length(levels.Temp) * length(levels.Sex)),
           xEnd      = rep(c(2,3,4,3,4,4) - 0.05, times = length(levels.Temp) * length(levels.Sex)),
           yLineRel  = rep((1 + c(0,1,3,0,2,0)) / 4.8, times = length(levels.Temp) * length(levels.Sex)),
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi)
  
}
```


## all plots month | temp * sex
```{r}
plots.MonthTempSex <- list()
tbl.plotResults.MonthTempSex <- data.frame(matrix(ncol = 6, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  dataMax <- listResults.MonthTempSex.plotPara[[m]]$dataMax
  yLo <- listResults.MonthTempSex.plotPara[[m]]$yLo
  yHi <- listResults.MonthTempSex.plotPara[[m]]$yHi
  
  plots.MonthTempSex[[m]] <- ggplot(data = listResults.MonthTempSex.means[[m]])
  
  # define layout parameters
  plots.MonthTempSex[[m]] <- plots.MonthTempSex[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
  facet_grid(rows = vars(sex), cols = vars(tempLevel), labeller = labeller(sex = labels.Sex))
  
  # create plot
  plots.MonthTempSex[[m]] <- plots.MonthTempSex[[m]] +
    geom_bar(mapping = aes(x = monthGroup,
                           y = emmean,
                           fill = monthGroup),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.MonthTempSex[[m]] <- plots.MonthTempSex[[m]] +
    geom_errorbar(mapping = aes(x = monthGroup,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.MonthTempSex[[m]] <- plots.MonthTempSex[[m]] +
    scale_x_discrete(guide = "axis_nested") # guide type comes from package "ggh4x"
  
  # add p-values
  plots.MonthTempSex[[m]] <- plots.MonthTempSex[[m]] +
    geom_segment(data = listResults.MonthTempSex.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.MonthTempSex.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.MonthTempSex[[m]] <- plots.MonthTempSex[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.MonthTempSex[[m]] <- plots.MonthTempSex[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.Temp,
       caption = paste("Model:", format(listFormula.full[[m]])),
       x = title.Month,
       y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.MonthTempSex[[m]] <- plots.MonthTempSex[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_MonthTempSex", ".svg")),
         width = 16, height = 16, units = "cm", dpi = 600, device = "svg",
         plot = plots.MonthTempSex[[m]], scale = 1)
  
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_MonthTempSex", ".png")),
         width = 16, height = 16, units = "cm", dpi = 600, device = "png",
         plot = plots.MonthTempSex[[m]], scale = 1)
  
  # results table
  tbl.plotResults.MonthTempSex <- rbind(tbl.plotResults.MonthTempSex, data.frame(rep(listOutcome[m], times = nrow(listResults.MonthTempSex.means[[1]])), listResults.MonthTempSex.means[[m]]$monthGroup, listResults.MonthTempSex.means[[m]]$tempLevel,  listResults.MonthTempSex.means[[m]]$sex, listResults.MonthTempSex.means[[m]]$emmean, listResults.MonthTempSex.means[[m]]$SE, listResults.MonthTempSex.means[[m]]$upper.CL - listResults.MonthTempSex.means[[m]]$emmean))
}

colnames(tbl.plotResults.MonthTempSex) <- c("targetVar", "monthGroup", "tempLevel", "sex", "mean", "SEM", "CI")
write.csv(tbl.plotResults.MonthTempSex, file.path(pathWD, paste0("tblPlotResults", "_MonthTempSex", ".csv")), row.names = FALSE)

plots.MonthTempSex

```


# temp | month * sex, results
```{r}
listResults.TempMonthSex <- purrr::map(listModel.final,
                                       ~ emmeans(.x,
                                                 pairwise ~ tempLevel | monthGroup * sex,
                                                 lmer.df = "satterthwaite",
                                                 adjust = NULL,
                                                 lmerTest.limit = 1e6)
)

listResults.TempMonthSex.means <- purrr::map(listResults.TempMonthSex,
                                             ~ .x$emmeans %>% as.data.frame()
)
listResults.TempMonthSex.means <- purrr::map(listResults.TempMonthSex.means,
                                             ~ within(.x, {
                                               lower.SEM <- emmean - SE
                                               upper.SEM <- emmean + SE})
)

listResults.TempMonthSex.contrasts <- purrr::map(listResults.TempMonthSex,
                                                 ~ .x$contrasts %>% as.data.frame()
)
```

### plot preparation
```{r}
yMinFrac <- c(0.4, 0.1, 0.1, 0.1, 0.1, 0.3, 0.1)
sigFrac <- 0.2

# adapt contrasts data set for plotting
listResults.TempMonthSex.plotPara <- list()
n = nrow(listResults.TempMonthSex.contrasts[[1]])
for (m in 1:length(listModel.final)) {

  dataMax <- listResults.TempMonthSex.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.TempMonthSex.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo     = yLo,
                                               yHi     = yHi)
  
  listResults.TempMonthSex.contrasts[[m]] <- listResults.TempMonthSex.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1,1,2) + 0.05, times = length(levels.Month) * length(levels.Sex)),
           xEnd     = rep(c(2,3,3) - 0.05, times = length(levels.Month) * length(levels.Sex)),
           yLineRel = rep((1 + c(0,1,0) ) / 2.8, times = length(levels.Month) * length(levels.Sex)),
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi)

}
```


## all plots temp | month * sex
```{r}
plots.TempMonthSex <- list()
tbl.plotResults.TempMonthSex <- data.frame(matrix(ncol = 6, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  dataMax <- listResults.TempMonthSex.plotPara[[m]]$dataMax
  yLo <- listResults.TempMonthSex.plotPara[[m]]$yLo
  yHi <- listResults.TempMonthSex.plotPara[[m]]$yHi
  
  plots.TempMonthSex[[m]] <- ggplot(data = listResults.TempMonthSex.means[[m]])
  
  # define layout parameters
  plots.TempMonthSex[[m]] <- plots.TempMonthSex[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
  facet_grid(rows = vars(sex), cols = vars(monthGroup), labeller = labeller(sex = labels.Sex))
  
  # create plot
  plots.TempMonthSex[[m]] <- plots.TempMonthSex[[m]] +
    geom_bar(mapping = aes(x = tempLevel,
                           y = emmean,
                           fill = tempLevel),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.TempMonthSex[[m]] <- plots.TempMonthSex[[m]] +
    geom_errorbar(mapping = aes(x = tempLevel,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  # plots.TempMonthSex[[m]] <- plots.TempMonthSex[[m]] +
  #   scale_x_discrete(guide = "axis_nested") # guide type comes from package "ggh4x"
  
  # add p-values
  plots.TempMonthSex[[m]] <- plots.TempMonthSex[[m]] +
    geom_segment(data = listResults.TempMonthSex.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.TempMonthSex.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.TempMonthSex[[m]] <- plots.TempMonthSex[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.TempMonthSex[[m]] <- plots.TempMonthSex[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.Month,
       caption = paste("Model:", format(listFormula.full[[m]])),
       x = title.Temp,
       y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.TempMonthSex[[m]] <- plots.TempMonthSex[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempMonthSex", ".svg")),
         width = 16, height = 16, units = "cm", dpi = 600, device = "svg",
         plot = plots.TempMonthSex[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempMonthSex", ".png")),
         width = 16, height = 16, units = "cm", dpi = 600, device = "png",
         plot = plots.TempMonthSex[[m]], scale = 1)
  
  # results table
  tbl.plotResults.TempMonthSex <- rbind(tbl.plotResults.TempMonthSex, data.frame(rep(listOutcome[m], times = nrow(listResults.TempMonthSex.means[[1]])), listResults.TempMonthSex.means[[m]]$tempLevel, listResults.TempMonthSex.means[[m]]$monthGroup,  listResults.TempMonthSex.means[[m]]$sex, listResults.TempMonthSex.means[[m]]$emmean, listResults.TempMonthSex.means[[m]]$SE, listResults.TempMonthSex.means[[m]]$upper.CL - listResults.TempMonthSex.means[[m]]$emmean))
}

colnames(tbl.plotResults.TempMonthSex) <- c("targetVar", "tempLevel", "monthGroup", "sex", "mean", "SEM", "CI")
write.csv(tbl.plotResults.TempMonthSex, file.path(pathWD, paste0("tblPlotResults", "_TempMonthSex", ".csv")), row.names = FALSE)

plots.TempMonthSex

```


# sex | month * temp, results
```{r}
listResults.SexMonthTemp <- purrr::map(listModel.final,
                                       ~ emmeans(.x,
                                                 pairwise ~ sex | monthGroup * tempLevel,
                                                 lmer.df = "satterthwaite",
                                                 adjust = NULL,
                                                 lmerTest.limit = 1e6)
)

listResults.SexMonthTemp.means <- purrr::map(listResults.SexMonthTemp,
                                             ~ .x$emmeans %>% as.data.frame()
)
listResults.SexMonthTemp.means <- purrr::map(listResults.SexMonthTemp.means,
                                             ~ within(.x, {
                                               lower.SEM <- emmean - SE
                                               upper.SEM <- emmean + SE})
)

listResults.SexMonthTemp.contrasts <- purrr::map(listResults.SexMonthTemp,
                                                 ~ .x$contrasts %>% as.data.frame()
)
```

### plot preparation
```{r}
yMinFrac <- c(0.4, 0.1, 0.1, 0.1, 0.1, 0.0, 0.0)
sigFrac <- 0.2

# adapt contrasts data set for plotting
listResults.SexMonthTemp.plotPara <- list()
n = nrow(listResults.SexMonthTemp.contrasts[[1]])
for (m in 1:length(listModel.final)) {

  dataMax <- listResults.SexMonthTemp.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.SexMonthTemp.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo     = yLo,
                                               yHi     = yHi)
  
  listResults.SexMonthTemp.contrasts[[m]] <- listResults.SexMonthTemp.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1) + 0.05, times = length(levels.Month) * length(levels.Temp)),
           xEnd     = rep(c(2) - 0.05, times = length(levels.Month) * length(levels.Temp)),
           yLineRel = rep((1 + c(0) ) / 1.8, times = length(levels.Month) * length(levels.Temp)),
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi)

}
```


## all plots sex | month * temp
```{r}
plots.SexMonthTemp <- list()
tbl.plotResults.SexMonthTemp <- data.frame(matrix(ncol = 6, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  dataMax <- listResults.SexMonthTemp.plotPara[[m]]$dataMax
  yLo <- listResults.SexMonthTemp.plotPara[[m]]$yLo
  yHi <- listResults.SexMonthTemp.plotPara[[m]]$yHi
  
  plots.SexMonthTemp[[m]] <- ggplot(data = listResults.SexMonthTemp.means[[m]])
  
  # define layout parameters
  plots.SexMonthTemp[[m]] <- plots.SexMonthTemp[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
  facet_grid(rows = vars(tempLevel), cols = vars(monthGroup))
  
  # create plot
  plots.SexMonthTemp[[m]] <- plots.SexMonthTemp[[m]] +
    geom_bar(mapping = aes(x = sex,
                           y = emmean,
                           fill = sex),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.SexMonthTemp[[m]] <- plots.SexMonthTemp[[m]] +
    geom_errorbar(mapping = aes(x = sex,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.SexMonthTemp[[m]] <- plots.SexMonthTemp[[m]] +
    scale_x_discrete(labels = labels.Sex)
  
  # add p-values
  plots.SexMonthTemp[[m]] <- plots.SexMonthTemp[[m]] +
    geom_segment(data = listResults.SexMonthTemp.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.SexMonthTemp.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.SexMonthTemp[[m]] <- plots.SexMonthTemp[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.SexMonthTemp[[m]] <- plots.SexMonthTemp[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.Month,
         caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Sex,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.SexMonthTemp[[m]] <- plots.SexMonthTemp[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_SexMonthTemp", ".svg")),
         width = 16, height = 16, units = "cm", dpi = 600, device = "svg",
         plot = plots.SexMonthTemp[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_SexMonthTemp", ".png")),
         width = 16, height = 16, units = "cm", dpi = 600, device = "png",
         plot = plots.SexMonthTemp[[m]], scale = 1)
  
  # results table
  tbl.plotResults.SexMonthTemp <- rbind(tbl.plotResults.SexMonthTemp, data.frame(rep(listOutcome[m], times = nrow(listResults.SexMonthTemp.means[[1]])), listResults.SexMonthTemp.means[[m]]$sex, listResults.SexMonthTemp.means[[m]]$monthGroup,  listResults.SexMonthTemp.means[[m]]$tempLevel, listResults.SexMonthTemp.means[[m]]$emmean, listResults.SexMonthTemp.means[[m]]$SE, listResults.SexMonthTemp.means[[m]]$upper.CL - listResults.SexMonthTemp.means[[m]]$emmean))
}

colnames(tbl.plotResults.SexMonthTemp) <- c("targetVar", "sex", "monthGroup", "tempLevel", "mean", "SEM", "CI")
write.csv(tbl.plotResults.SexMonthTemp, file.path(pathWD, paste0("tblPlotResults", "_SexMonthTemp", ".csv")), row.names = FALSE)

plots.SexMonthTemp

```


# temp @ capture
## model
```{r}
# formula
formula.tempAtCapture.full <- as.formula(paste("tempAtCapture", "~ monthGroup + (1|replicate)"))

# calculate model
model.tempAtCapture.full <- lmer(formula = formula.tempAtCapture.full,
                                 data = combTable,
                                 REML = TRUE,
                                 na.action = na.exclude)

# step function from lmerTest package
# remove interactions which are not significant
# but keep always the main effect
model.tempAtCapture.step <- lmerTest::step(model.tempAtCapture.full,
                                           ddf = "Satterthwaite",
                                           alpha.fixed = 0.05,
                                           alpha.random = 0.1,
                                           reduce.fixed = TRUE,
                                           reduce.random = FALSE,
                                           keep = c("monthGroup")
                             )

# extract final model
model.tempAtCapture.final <- get_model(model.tempAtCapture.step)
```


## fitted values
```{r}
results.tempAtCapture.fittedValues <- data.frame(replicate = combTable$replicate,
                                                 monthGroup = combTable$monthGroup,
                                                 animalType = combTable$animalType,
                                                 sex = combTable$sex,
                                                 tempLevel = combTable$tempLevel,
                                                 fitted = fitted.values(model.tempAtCapture.final))

# remove NA
results.tempAtCapture.fittedValues <- na.omit(results.tempAtCapture.fittedValues)

# remove all duplicates
results.tempAtCapture.fittedValues <-  results.tempAtCapture.fittedValues[!duplicated(results.tempAtCapture.fittedValues[c('replicate','animalType','tempLevel')]), ]
```


### mean over tempLevel
```{r}
library("dplyr")

results.tempAtCapture.fittedValues.mean <- results.tempAtCapture.fittedValues %>%
  group_by(replicate,monthGroup,animalType,sex) %>% 
  summarise_at(vars(fitted), list(fitted.mean = mean))
```


## results
```{r}
results.tempAtCapture <- emmeans(model.tempAtCapture.final,
                                 pairwise ~ monthGroup,
                                 lmer.df = "satterthwaite",
                                 adjust = NULL,
                                 lmerTest.limit = 1e6)

results.tempAtCapture.means <- results.tempAtCapture$emmeans %>% as.data.frame()
results.tempAtCapture.means <- within(results.tempAtCapture.means, {
                                       lower.SEM <- emmean - SE
                                       upper.SEM <- emmean + SE})

results.tempAtCapture.contrasts <- results.tempAtCapture$contrasts %>% as.data.frame()

# adjustment for multiple comparisons: Tukey test
# results.tempAtCapture.contrasts <- within(results.tempAtCapture.contrasts,
#                                           p.value.adj <- 1 - ptukey(results.tempAtCapture.contrasts$t.ratio * sqrt(2),
#                                                                     nmeans = nrow(results.tempAtCapture.contrasts),
#                                                                     df = results.tempAtCapture.contrasts$df))
# adjustment for multiple comparisons: Benjamini Hochberg
results.tempAtCapture.contrasts <- within(results.tempAtCapture.contrasts,
                                          p.value.adj <- p.adjust(results.tempAtCapture.contrasts$p.value,
                                                                  method = "BH")
)

rm(results.tempAtCapture)
```

## plot preparation
```{r}
# adapt contrasts data set for plotting
results.tempAtCapture.contrasts <- results.tempAtCapture.contrasts %>%
  mutate(sig.idx = case_when(p.value.adj < 0.001 ~ "***",
                             p.value.adj < 0.01  ~ "**",
                             p.value.adj < 0.05  ~ "*",
                             TRUE ~ "")) %>%
  mutate(p.valuef = sprintf("%1.4f", p.value.adj)) %>%
  mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                              TRUE ~ as.character(p.valuef))) %>%
  mutate(xStart   = c(1,1,1,2,2,3) + 0.05,
         xEnd     = c(2,3,4,3,4,4) - 0.05,
         xText  = 0.5*(xStart+xEnd))
```



## plot LMM (point) v2
```{r fig.height=8/2.54, fig.width=9.3/2.54}
plots.TempAtCapture.limits <- data.frame(min = c(0),
                                          max = c(36))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  results.tempAtCapture.contrasts <- results.tempAtCapture.contrasts %>%
    mutate(yStart = rep((0.69 + 0.08*c(0,1,3,0,2,0)) * (plots.TempAtCapture.limits$max - plots.TempAtCapture.limits$min) + plots.TempAtCapture.limits$min, times = 1),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.TempAtCapture.limits$max - plots.TempAtCapture.limits$min))
}

tbl.plotResults.tempAtCapture <- data.frame(matrix(ncol = 4, nrow = 0))

plot.TempAtCapture <- ggplot(data = results.tempAtCapture.means)

# define layout parameters
plot.TempAtCapture <- plot.TempAtCapture +
  theme_bw() +
  theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
        axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
        axis.title.y = element_text(size = 15, margin = margin(r=10)),
        axis.text.y = element_text(size = 10, face = "bold"),
        panel.border = element_rect(size = 1.5, colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
        plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
        plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
        legend.position = "none")

# create point plot: replicate means
plot.TempAtCapture <- plot.TempAtCapture +
  geom_beeswarm(data = results.tempAtCapture.fittedValues.mean,
                mapping = aes(x = monthGroup,
                              y = fitted.mean),
                colour = "#C0C0C0",
                size = 2,
                cex = 1.6)
  
# create plot
plot.TempAtCapture <- plot.TempAtCapture +
  geom_point(mapping = aes(x = monthGroup,
                           y = emmean),
             stat = "identity", # important since we provide mean values already
             size = 3, # line width of bar borders
             colour = "#000000") # bar border color

# add error bars
plot.TempAtCapture <- plot.TempAtCapture +
  geom_errorbar(mapping = aes(x = monthGroup,
                              ymin = lower.SEM,
                              ymax = upper.SEM),
                size=0.75,
                width=0.25,
                colour = "#000000",
                position = position_dodge(0.1))

# add discrete x-axis labels
plot.TempAtCapture <- plot.TempAtCapture +
  scale_x_discrete(labels = labels.Month.Short)

# add p-values
plot.TempAtCapture <- plot.TempAtCapture +
  geom_segment(data = results.tempAtCapture.contrasts,
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = results.tempAtCapture.contrasts,
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)

# add nicer fill color
plot.TempAtCapture <- plot.TempAtCapture +
  scale_fill_brewer(palette = "Blues")

# add further layout parameter
plot.TempAtCapture <- plot.TempAtCapture +
  scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                     limits = c(0, plots.TempAtCapture.limits$max),
                     expand = expansion(mult = c(0.025, 0.025))) +
  labs(caption = paste("Model:", format(formula.tempAtCapture.full)),
       x = title.Month,
       y = title.Temp.Vicinity)

# change coordinates = "zoom in"
plot.TempAtCapture <- plot.TempAtCapture +
  coord_cartesian(ylim = c(plots.TempAtCapture.limits$min, plots.TempAtCapture.limits$max))

# save plot
ggsave(filename = file.path(pathWD, paste0("TempAtCapture", "_LMM", ".svg")),
       width = 9.3, height = 8, units = "cm", dpi = 600, device = "svg",
       plot = plot.TempAtCapture, scale = 1)

ggsave(filename = file.path(pathWD, paste0("TempAtCapture", "_LMM", ".png")),
       width = 9.3, height = 8, units = "cm", dpi = 600, device = "png",
       plot = plot.TempAtCapture, scale = 1)

# results table
  tbl.plotResults.tempAtCapture <- data.frame(rep(title.Temp.Vicinity, times = nrow(results.tempAtCapture.means)), results.tempAtCapture.means$monthGroup, results.tempAtCapture.means$emmean, results.tempAtCapture.means$SE, results.tempAtCapture.means$upper.CL - results.tempAtCapture.means$emmean)

colnames(tbl.plotResults.tempAtCapture) <- c("targetVar", "tempAtCapture", "mean", "SEM", "CI")
write.csv(tbl.plotResults.tempAtCapture, file.path(pathWD, paste0("tblPlotResults", "_tempAtCapture", ".csv")), row.names = FALSE)

plot.TempAtCapture

```

### save @ another size
```{r fig.height=10/2.54, fig.width=12/2.54}
# save plot
ggsave(filename = file.path(pathWD, paste0("TempAtCapture", "_LMM", "_big", ".svg")),
       width = 12, height = 10, units = "cm", dpi = 600, device = "svg",
       plot = plot.TempAtCapture, scale = 1)

ggsave(filename = file.path(pathWD, paste0("TempAtCapture", "_LMM", "_big", ".png")),
       width = 12, height = 10, units = "cm", dpi = 600, device = "png",
       plot = plot.TempAtCapture, scale = 1)

plot.TempAtCapture

```


## boxplot
```{r}
table.tempAtCapture <- combTable[ , c("monthGroup", "tempAtCapture", "replicate")]

# remove duplicates (individual cell measurements)
table.tempAtCapture <- table.tempAtCapture  %>%
  distinct(.keep_all = TRUE)


factorYMax <- 1.0
boxplot.TempAtCapture <- ggplot(data = table.tempAtCapture, aes(x = monthGroup, y = tempAtCapture, fill = monthGroup)) +
       stat_boxplot(geom = "errorbar", # Boxplot with error bars 
                    width = 0.2) +
       geom_boxplot(outlier.colour = "red")
# define layout parameters
boxplot.TempAtCapture <- boxplot.TempAtCapture +
  theme_bw() +
  theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
        axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
        axis.title.y = element_text(size = 15, margin = margin(r=10)),
        axis.text.y = element_text(size = 10, face = "bold"),
        panel.border = element_rect(size = 1.5, colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
        plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
        plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
        legend.position = "none")

# add nicer fill color
boxplot.TempAtCapture <- boxplot.TempAtCapture +
  scale_fill_brewer(palette = "Blues")

# add further layout parameter
ymax <- results.tempAtCapture.means %>% select(upper.CL) %>% max(., na.rm = TRUE)
boxplot.TempAtCapture <- boxplot.TempAtCapture +
  scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                     limits = c(0, factorYMax*ymax),
                     expand = expansion(mult = c(0.025, 0.025))) +
  labs(caption = "Boxplot",
       x = title.Month,
       y = title.Temp.Vicinity)

# change coordinates = "zoom in"
boxplot.TempAtCapture <- boxplot.TempAtCapture +
  coord_cartesian(ylim = c(ymax*.0, factorYMax*ymax))

# save plot
ggsave(filename = file.path(pathWD, paste0("TempAtCapture", "_Box", ".svg")),
       width = 10, height = 8, units = "cm", dpi = 600, device = "svg",
       plot = boxplot.TempAtCapture, scale = 1)

ggsave(filename = file.path(pathWD, paste0("TempAtCapture", "_Box", ".png")),
       width = 10, height = 8, units = "cm", dpi = 600, device = "png",
       plot = boxplot.TempAtCapture, scale = 1)

# save plot for presentation
# ggsave(filename = file.path(pathWD, paste0("TempAtCapture", "_Box_pres", ".svg")),
#        width = 11, height = 8, units = "cm", dpi = 600, device = "svg",
#        plot = boxplot.TempAtCapture, scale = 1)
# 
# ggsave(filename = file.path(pathWD, paste0("TempAtCapture", "_Box_pres", ".png")),
#        width = 11, height = 8, units = "cm", dpi = 600, device = "png",
#        plot = boxplot.TempAtCapture, scale = 1)

boxplot.TempAtCapture

```


# combined plot: month & temp @ capture
## prepare data (combine)
```{r}
listResults.Month_and_TempAtCapture <- list()

for (m in 1:length(listModel.final)) {
  listResults.Month_and_TempAtCapture[[m]] <- cbind(listResults.Month.means[[m]], 
                                                    data.frame(tempAtCptr_mean = results.tempAtCapture.means$emmean,
                                                               tempAtCptr_SE = results.tempAtCapture.means$SE,
                                                               tempAtCptr_df = results.tempAtCapture.means$df,
                                                               tempAtCptr_lower.CL = results.tempAtCapture.means$lower.CL,
                                                               tempAtCptr_upper.CL = results.tempAtCapture.means$upper.CL,
                                                               tempAtCptr_lower.SEM = results.tempAtCapture.means$lower.SEM,
                                                               tempAtCptr_upper.SEM = results.tempAtCapture.means$upper.SEM))
  y1Lo <- plots.Month.limits$min[m]
  y1Hi <- plots.Month.limits$max[m]
  
  y2Lo <- plots.TempAtCapture.limits$min
  y2Hi <- plots.TempAtCapture.limits$max
  
  # scaling between axes
  b <- (y1Hi - y1Lo) / (y2Hi - y2Lo)
  a <- y1Lo - b*y2Lo
  
  listResults.Month_and_TempAtCapture[[m]] <- cbind(listResults.Month_and_TempAtCapture[[m]],
                                                    data.frame(scaleOffset = a,
                                                               scaleSlope = b,
                                                               tempAtCptr_scaledMean = a + b*listResults.Month_and_TempAtCapture[[m]]$tempAtCptr_mean,
                                                               tempAtCptr_scaledLower.CL = a + b*listResults.Month_and_TempAtCapture[[m]]$tempAtCptr_lower.CL,
                                                               tempAtCptr_scaledUpper.CL = a + b*listResults.Month_and_TempAtCapture[[m]]$tempAtCptr_upper.CL,
                                                               tempAtCptr_scaledLower.SEM = a + b*listResults.Month_and_TempAtCapture[[m]]$tempAtCptr_lower.SEM,
                                                               tempAtCptr_scaledUpper.SEM = a + b*listResults.Month_and_TempAtCapture[[m]]$tempAtCptr_upper.SEM))
}
rm(y1Lo, y1Hi, y2Lo, y2Hi, a, b)
```


## all plots
```{r}
plots.Month_and_TempAtCapture <- list()

for (m in 1:length(listModel.final)) {
  
  y1Lo <- plots.Month.limits$min[m]
  y1Hi <- plots.Month.limits$max[m]
  
  scaleOffset <- listResults.Month_and_TempAtCapture[[m]]$scaleOffset[1]
  scaleSlope <- listResults.Month_and_TempAtCapture[[m]]$scaleSlope[1]
  
  plots.Month_and_TempAtCapture[[m]] <- ggplot(data = listResults.Month_and_TempAtCapture[[m]], aes(x = monthGroup))
  
  # define layout parameters
  plots.Month_and_TempAtCapture[[m]] <- plots.Month_and_TempAtCapture[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10), color = "#1f77bd"),
          axis.text.y = element_text(size = 10, face = "bold"),
          axis.title.y.right = element_text(size = 15, face = "plain", margin = margin(l=10), color = "#c21618"),
          axis.text.y.right = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")
  
  # create plot
  plots.Month_and_TempAtCapture[[m]] <- plots.Month_and_TempAtCapture[[m]] +
    geom_point(mapping = aes(y = emmean),
             stat = "identity", # important since we provide mean values already
             size = 3, # line width of bar borders
             colour = "#1f77bd") # bar border color
  
  # add error bars
  plots.Month_and_TempAtCapture[[m]] <- plots.Month_and_TempAtCapture[[m]] +
    geom_errorbar(mapping = aes(ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = "#1f77bd",
                  position = position_dodge(0.1))
  
  # add nested x-axis
  # plots.Month_and_TempAtCapture[[m]] <- plots.Month_and_TempAtCapture[[m]] +
  #   scale_x_discrete(guide = "axis_nested") # guide type comes from package "ggh4x"
  
  # add p-values
  plots.Month_and_TempAtCapture[[m]] <- plots.Month_and_TempAtCapture[[m]] +
    geom_segment(data = listResults.Month.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd),
                 size = 0.5, colour="#1f77bd", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.Month.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7,
              colour = "#1f77bd")
  
  # temperature @ capture as second axis
  plots.Month_and_TempAtCapture[[m]] <- plots.Month_and_TempAtCapture[[m]] +
    geom_point(mapping = aes(y = tempAtCptr_scaledMean),
             stat = "identity", # important since we provide mean values already
             size = 3, # line width of bar borders
             colour = "#c21618") # bar border color
  
  # add error bars
  plots.Month_and_TempAtCapture[[m]] <- plots.Month_and_TempAtCapture[[m]] +
    geom_errorbar(mapping = aes(ymin = tempAtCptr_scaledLower.SEM,
                                ymax = tempAtCptr_scaledUpper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = "#c21618",
                  position = position_dodge(0.1))
  
  # add further layout parameter
  plots.Month_and_TempAtCapture[[m]] <- plots.Month_and_TempAtCapture[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(y1Lo, y1Hi),
                       expand = expansion(mult = c(0.025, 0.025)),
                       sec.axis = sec_axis(~ (. - scaleOffset) / scaleSlope, name = title.Temp.Vicinity)) +
    labs(caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Month,
         y = listLabelOutcome[m])
  
  plots.Month_and_TempAtCapture[[m]] %>% print()
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Month_and_TempAtCapture", ".svg")),
         width = 10, height = 8, units = "cm", dpi = 600, device = "svg",
         plot = plots.Month_and_TempAtCapture[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Month_and_TempAtCapture", ".png")),
         width = 10, height = 8, units = "cm", dpi = 600, device = "png",
         plot = plots.Month_and_TempAtCapture[[m]], scale = 1)
}

```


# combined plot: temp | month & temp @ capture
## prepare data (combine)
```{r}
listResults.TempMonth_and_TempAtCapture <- list()

for (m in 1:length(listModel.final)) {
  listResults.TempMonth_and_TempAtCapture[[m]] <- cbind(listResults.TempMonth.means[[m]], 
                                                        data.frame(tempAtCptr_mean = rep(results.tempAtCapture.means$emmean, each = length(levels.Temp)) * rep(c(NaN,1,NaN), times = length(levels.Month)),
                                                                   tempAtCptr_SE = rep(results.tempAtCapture.means$SE, each = length(levels.Temp)) * rep(c(NaN,1,NaN), times = length(levels.Month)),
                                                                   tempAtCptr_df = rep(results.tempAtCapture.means$df, each = length(levels.Temp)) * rep(c(NaN,1,NaN), times = length(levels.Month)),
                                                                   tempAtCptr_lower.CL = rep(results.tempAtCapture.means$lower.CL, each = length(levels.Temp)) * rep(c(NaN,1,NaN), times = length(levels.Month)),
                                                                   tempAtCptr_upper.CL = rep(results.tempAtCapture.means$upper.CL, each = length(levels.Temp)) * rep(c(NaN,1,NaN), times = length(levels.Month)),
                                                                   tempAtCptr_lower.SEM = rep(results.tempAtCapture.means$lower.SEM, each = length(levels.Temp)) * rep(c(NaN,1,NaN), times = length(levels.Month)),
                                                                   tempAtCptr_upper.SEM = rep(results.tempAtCapture.means$upper.SEM, each = length(levels.Temp)) * rep(c(NaN,1,NaN), times = length(levels.Month))
                                                                   )
                                                        )
  
  y1Lo <- listResults.TempMonth.plotPara[[m]]$yLo
  y1Hi <- listResults.TempMonth.plotPara[[m]]$yHi
  
  y2Lo <- plots.TempAtCapture.limits$min
  y2Hi <- plots.TempAtCapture.limits$max
  
  # scaling between axes
  b <- (y1Hi - y1Lo) / (y2Hi - y2Lo)
  a <- y1Lo - b*y2Lo
  
  listResults.TempMonth_and_TempAtCapture[[m]] <- cbind(listResults.TempMonth_and_TempAtCapture[[m]],
                                                        data.frame(scaleOffset = a,
                                                                   scaleSlope = b,
                                                                   tempAtCptr_scaledMean = a + b*listResults.TempMonth_and_TempAtCapture[[m]]$tempAtCptr_mean,
                                                                   tempAtCptr_scaledLower.CL = a + b*listResults.TempMonth_and_TempAtCapture[[m]]$tempAtCptr_lower.CL,
                                                                   tempAtCptr_scaledUpper.CL = a + b*listResults.TempMonth_and_TempAtCapture[[m]]$tempAtCptr_upper.CL,
                                                                   tempAtCptr_scaledLower.SEM = a + b*listResults.TempMonth_and_TempAtCapture[[m]]$tempAtCptr_lower.SEM,
                                                                   tempAtCptr_scaledUpper.SEM = a + b*listResults.TempMonth_and_TempAtCapture[[m]]$tempAtCptr_upper.SEM))
}
rm(y1Lo, y1Hi, y2Lo, y2Hi, a, b)
```


## all plots
```{r}
plots.TempMonth_and_TempAtCapture <- list()

for (m in 1:length(listModel.final)) {
  
  y1Lo <- listResults.TempMonth.plotPara[[m]]$yLo
  y1Hi <- listResults.TempMonth.plotPara[[m]]$yHi
  
  scaleOffset <- listResults.TempMonth_and_TempAtCapture[[m]]$scaleOffset[1]
  scaleSlope <- listResults.TempMonth_and_TempAtCapture[[m]]$scaleSlope[1]
  
  plots.TempMonth_and_TempAtCapture[[m]] <- ggplot(data = listResults.TempMonth_and_TempAtCapture[[m]], aes(x = interaction(tempLevel, monthGroup)))
  
  # define layout parameters
  plots.TempMonth_and_TempAtCapture[[m]] <- plots.TempMonth_and_TempAtCapture[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10), color = "#1f77bd"),
          axis.text.y = element_text(size = 10, face = "bold"),
          axis.title.y.right = element_text(size = 15, face = "plain", margin = margin(l=10), color = "#c21618"),
          axis.text.y.right = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")
  
  # create plot
  plots.TempMonth_and_TempAtCapture[[m]] <- plots.TempMonth_and_TempAtCapture[[m]] +
    geom_point(mapping = aes(y = emmean),
             stat = "identity", # important since we provide mean values already
             size = 3, # line width of bar borders
             colour = "#1f77bd") # bar border color
  
  # add error bars
  plots.TempMonth_and_TempAtCapture[[m]] <- plots.TempMonth_and_TempAtCapture[[m]] +
    geom_errorbar(mapping = aes(ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = "#1f77bd",
                  position = position_dodge(0.1))
  
  # add nested x-axis
  plots.TempMonth_and_TempAtCapture[[m]] <- plots.TempMonth_and_TempAtCapture[[m]] +
    scale_x_discrete(guide = "axis_nested") # guide type comes from package "ggh4x"
  
  # add p-values
  plots.TempMonth_and_TempAtCapture[[m]] <- plots.TempMonth_and_TempAtCapture[[m]] +
    geom_segment(data = listResults.TempMonth.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd),
                 size = 0.5, colour="#1f77bd", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.TempMonth.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7,
              colour = "#1f77bd")
  
  # temperature @ capture as second axis
  plots.TempMonth_and_TempAtCapture[[m]] <- plots.TempMonth_and_TempAtCapture[[m]] +
    geom_point(mapping = aes(y = tempAtCptr_scaledMean),
             stat = "identity", # important since we provide mean values already
             size = 3, # line width of bar borders
             colour = "#c21618") # bar border color

  # add error bars
  plots.TempMonth_and_TempAtCapture[[m]] <- plots.TempMonth_and_TempAtCapture[[m]] +
    geom_errorbar(mapping = aes(ymin = tempAtCptr_scaledLower.SEM,
                                ymax = tempAtCptr_scaledUpper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = "#c21618",
                  position = position_dodge(0.1))

  # add further layout parameter
  plots.TempMonth_and_TempAtCapture[[m]] <- plots.TempMonth_and_TempAtCapture[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(y1Lo, y1Hi),
                       expand = expansion(mult = c(0.025, 0.025)),
                       sec.axis = sec_axis(~ (. - scaleOffset) / scaleSlope, name = title.Temp.Vicinity)) +
    labs(caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Month,
         y = listLabelOutcome[m])

  plots.TempMonth_and_TempAtCapture[[m]] %>% print()

  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempMonth_and_TempAtCapture", ".svg")),
         width = 16, height = 8, units = "cm", dpi = 600, device = "svg",
         plot = plots.TempMonth_and_TempAtCapture[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempMonth_and_TempAtCapture", ".png")),
         width = 16, height = 8, units = "cm", dpi = 600, device = "png",
         plot = plots.TempMonth_and_TempAtCapture[[m]], scale = 1)
}

```


# effect size, reference levels
```{r}
refLevel.effectSize.month <- "4-5"

# labels for outcome variables: effect size
listLabelEffectSize <- c("Relative cell size",
                         expression("Relative " ~ tau["inlet"]),
                         expression("Relative " ~ tau["channel"]),
                         expression("Relative " ~ hat(d)["inlet"]),
                         expression("Relative " ~ hat(d)["channel"]),
                         "Relative Young's modulus",
                         "Relative viscosity")

```


# effect size: month, results
```{r}
listEffectSizes.Month.means <- list()
listEffectSizes.Month.contrasts <- list()

for (m in 1:length(listModel.final)) {
  # copy contrast data
  listEffectSizes.Month.contrasts[[m]] <- data.frame(contrast = listResults.Month.contrasts[[m]]$contrast,
                                                     df = listResults.Month.contrasts[[m]]$df,
                                                     t.ratio = listResults.Month.contrasts[[m]]$t.ratio,
                                                     p.value = listResults.Month.contrasts[[m]]$p.value,
                                                     p.value.adj = listResults.Month.contrasts[[m]]$p.value.adj,
                                                     sig.idx = listResults.Month.contrasts[[m]]$sig.idx,
                                                     p.valuef = listResults.Month.contrasts[[m]]$p.valuef)
  
  # copy column of monthGroup to results list for effect size
  listEffectSizes.Month.means[[m]] <- data.frame(monthGroup = listResults.Month.means[[m]]$monthGroup)
  # create columns for effect size and errors
  listEffectSizes.Month.means[[m]]$mean <- NaN
  listEffectSizes.Month.means[[m]]$SE <- NaN
  listEffectSizes.Month.means[[m]]$lower.CL <- NaN
  listEffectSizes.Month.means[[m]]$upper.CL <- NaN
  listEffectSizes.Month.means[[m]]$lower.SEM <- NaN
  listEffectSizes.Month.means[[m]]$upper.SEM <- NaN
  
  for (idx in 1:nrow(listEffectSizes.Month.means[[m]])) {
    if (listEffectSizes.Month.means[[m]]$monthGroup[idx] == refLevel.effectSize.month) {
      # reference
      listEffectSizes.Month.means[[m]]$mean[idx] <- 1
      listEffectSizes.Month.means[[m]]$SE[idx] <- listResults.Month.means[[m]]$SE[idx] / listResults.Month.means[[m]]$emmean[idx]
      listEffectSizes.Month.means[[m]]$lower.CL[idx] <- listResults.Month.means[[m]]$lower.CL[idx] / listResults.Month.means[[m]]$emmean[idx]
      listEffectSizes.Month.means[[m]]$upper.CL[idx] <- listResults.Month.means[[m]]$upper.CL[idx] / listResults.Month.means[[m]]$emmean[idx]
      listEffectSizes.Month.means[[m]]$lower.SEM[idx] <- listResults.Month.means[[m]]$lower.SEM[idx] / listResults.Month.means[[m]]$emmean[idx]
      listEffectSizes.Month.means[[m]]$upper.SEM[idx] <- listResults.Month.means[[m]]$upper.SEM[idx] / listResults.Month.means[[m]]$emmean[idx]
    } else {
      # other level: calculate uncertainty by linear error term analysis
      # find reference data
      refMask <- (listEffectSizes.Month.means[[m]]$monthGroup == refLevel.effectSize.month)
      # means of reference level and recent data row
      refLevel <- listResults.Month.means[[m]]$emmean[refMask]
      dataLevel <- listResults.Month.means[[m]]$emmean[idx]
      
      # error = confidence level
      refDelta <- listResults.Month.means[[m]]$upper.CL[refMask] - listResults.Month.means[[m]]$emmean[refMask]
      dataDelta <- listResults.Month.means[[m]]$upper.CL[idx] - listResults.Month.means[[m]]$emmean[idx]
      # effect size
      effectSize <- dataLevel / refLevel
      effectSizeDelta <- effectSize * sqrt( (dataDelta / dataLevel)^2 + (refDelta / refLevel)^2 )
      # assign results
      listEffectSizes.Month.means[[m]]$mean[idx] <- effectSize
      listEffectSizes.Month.means[[m]]$lower.CL[idx] <- effectSize - effectSizeDelta
      listEffectSizes.Month.means[[m]]$upper.CL[idx] <- effectSize + effectSizeDelta
      
      # error = SEM
      refDelta <- listResults.Month.means[[m]]$SE[refMask]
      dataDelta <- listResults.Month.means[[m]]$SE[idx]
      # effect size
      effectSizeDelta <- effectSize * sqrt( (dataDelta / dataLevel)^2 + (refDelta / refLevel)^2 )
      # assign results
      listEffectSizes.Month.means[[m]]$SE[idx] <- effectSizeDelta
      listEffectSizes.Month.means[[m]]$lower.SEM[idx] <- effectSize - effectSizeDelta
      listEffectSizes.Month.means[[m]]$upper.SEM[idx] <- effectSize + effectSizeDelta
    }
  }
}
rm(refMask, refLevel, dataLevel, refDelta, dataDelta, effectSize, effectSizeDelta)
```


## effect size: fitted values
```{r}
listEffectSizes.Month.fittedValues <- listResults.fittedValues.mean

for (m in 1:length(listModel.final)) {
  listEffectSizes.Month.fittedValues[[m]] <- listEffectSizes.Month.fittedValues[[m]] %>%
    mutate(fitted.normalized = NA)
  
  for (idx in 1:nrow(listEffectSizes.Month.fittedValues[[m]])) {
    refMask <- (listResults.Month.means[[m]]$monthGroup == refLevel.effectSize.month)
    # means of reference level and recent data row
    refLevel <- listResults.Month.means[[m]]$emmean[refMask]
      
    listEffectSizes.Month.fittedValues[[m]]$fitted.normalized[idx] <- listEffectSizes.Month.fittedValues[[m]]$fitted.mean[idx] / refLevel
  }
}

```


## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listEffectSizes.Month.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listEffectSizes.Month.contrasts[[m]] <- listEffectSizes.Month.contrasts[[m]] %>%
    mutate(xStart   = c(1,1,1,2,2,3) + 0.05,
           xEnd     = c(2,3,4,3,4,4) - 0.05,
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots effect size: month (point) v2
```{r fig.height=6/2.54, fig.width=9.3/2.54}
plots.effectSize.Month.limits <- data.frame(min = c(.875, 0.5, 0.5, 0.5, 0.5, 0.9, 0.6),
                                            max = c(1.125, 2.3, 2.3, 1.7, 1.7, 1.05, 1.3))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listEffectSizes.Month.contrasts[[m]] <- listEffectSizes.Month.contrasts[[m]] %>%
    mutate(yStart = rep((0.69 + 0.08*c(0,1,3,0,2,0)) * (plots.effectSize.Month.limits$max[m] - plots.effectSize.Month.limits$min[m]) + plots.effectSize.Month.limits$min[m], times = 1),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.effectSize.Month.limits$max[m] - plots.effectSize.Month.limits$min[m]))
}

plots.effectSize.Month <- list()
tbl.plotResults.effectSize.Month <- data.frame(matrix(ncol = 4, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.effectSize.Month[[m]] <- ggplot(data = listEffectSizes.Month.means[[m]])
  
  # define layout parameters
  plots.effectSize.Month[[m]] <- plots.effectSize.Month[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")
  
  # create point plot: replicate means
  plots.effectSize.Month[[m]] <- plots.effectSize.Month[[m]] +
    geom_beeswarm(data = listEffectSizes.Month.fittedValues[[m]],
                  mapping = aes(x = monthGroup,
                                y = fitted.normalized),
                  colour = "#C0C0C0",
                  size = 2,
                  cex = 1.6)
  # are some means out of viewing window?
  if( (plots.effectSize.Month.limits$min[m] > min(listEffectSizes.Month.fittedValues[[m]]$fitted.normalized)) | (plots.effectSize.Month.limits$max[m] < max(listEffectSizes.Month.fittedValues[[m]]$fitted.normalized)) ) {
    print(paste("Replicate means out of view window:",listLabelEffectSize[m]))
  }
  
  # create plot
  plots.effectSize.Month[[m]] <- plots.effectSize.Month[[m]] +
    geom_point(mapping = aes(x = monthGroup,
                             y = mean),
             stat = "identity", # important since we provide mean values already
             size = 3, # line width of bar borders
             colour = "#1f77bd") # bar border color
  
  # add error bars
  plots.effectSize.Month[[m]] <- plots.effectSize.Month[[m]] +
    geom_errorbar(mapping = aes(x = monthGroup,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = "#1f77bd",
                  position = position_dodge(0.1))
  
  # add discrete x-axis labels
  plots.effectSize.Month[[m]] <- plots.effectSize.Month[[m]] +
    scale_x_discrete(labels = labels.Month.Short)

  # add p-values
  # plots.effectSize.Month[[m]] <- plots.effectSize.Month[[m]] +
  #   geom_segment(data = listEffectSizes.Month.contrasts[[m]],
  #                mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
  #                size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
  #   geom_text(data = listEffectSizes.Month.contrasts[[m]],
  #             mapping = aes(x = xText, y = yText, label = sig.idx),
  #             size = 7)
  
  # add nicer fill color
  plots.effectSize.Month[[m]] <- plots.effectSize.Month[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.effectSize.Month[[m]] <- plots.effectSize.Month[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.effectSize.Month.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Month,
         y = listLabelEffectSize[m])
  
  # change coordinates = "zoom in"
  plots.effectSize.Month[[m]] <- plots.effectSize.Month[[m]] +
    coord_cartesian(ylim = c(plots.effectSize.Month.limits$min[m], plots.effectSize.Month.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_effectSize", "_Month", ".svg")),
         width = 9.3, height = 6, units = "cm", dpi = 600, device = "svg",
         plot = plots.effectSize.Month[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_effectSize", "_Month", ".png")),
         width = 9.3, height = 6, units = "cm", dpi = 600, device = "png",
         plot = plots.effectSize.Month[[m]], scale = 1)
  
  # results table
  tbl.plotResults.effectSize.Month <- rbind(tbl.plotResults.effectSize.Month, data.frame(rep(listOutcome[m], times = nrow(listEffectSizes.Month.means[[1]])), listEffectSizes.Month.means[[m]]$monthGroup, listEffectSizes.Month.means[[m]]$mean, listEffectSizes.Month.means[[m]]$SE, listEffectSizes.Month.means[[m]]$upper.CL - listEffectSizes.Month.means[[m]]$mean))
}

colnames(tbl.plotResults.effectSize.Month) <- c("targetVar", "monthGroup", "mean", "SEM", "CI")
write.csv(tbl.plotResults.effectSize.Month, file.path(pathWD, paste0("tblPlotResults", "_effectSize", "_Month", ".csv")), row.names = FALSE)

plots.effectSize.Month

```


# effect size: month | temp, results
```{r}
listEffectSizes.MonthTemp.means <- list()
listEffectSizes.MonthTemp.contrasts <- list()

for (m in 1:length(listModel.final)) {
  # copy contrast data
  listEffectSizes.MonthTemp.contrasts[[m]] <- data.frame(contrast = listResults.MonthTemp.contrasts[[m]]$contrast,
                                                     tempLevel = listResults.MonthTemp.contrasts[[m]]$tempLevel,
                                                     df = listResults.MonthTemp.contrasts[[m]]$df,
                                                     t.ratio = listResults.MonthTemp.contrasts[[m]]$t.ratio,
                                                     p.value = listResults.MonthTemp.contrasts[[m]]$p.value,
                                                     p.value.adj = listResults.MonthTemp.contrasts[[m]]$p.value.adj,
                                                     sig.idx = listResults.MonthTemp.contrasts[[m]]$sig.idx,
                                                     p.valuef = listResults.MonthTemp.contrasts[[m]]$p.valuef)
  
  # copy column of monthGroup to results list for effect size
  listEffectSizes.MonthTemp.means[[m]] <- data.frame(monthGroup = listResults.MonthTemp.means[[m]]$monthGroup,
                                                     tempLevel = listResults.MonthTemp.means[[m]]$tempLevel)
  # create columns for effect size and errors
  listEffectSizes.MonthTemp.means[[m]]$mean <- NaN
  listEffectSizes.MonthTemp.means[[m]]$SE <- NaN
  listEffectSizes.MonthTemp.means[[m]]$lower.CL <- NaN
  listEffectSizes.MonthTemp.means[[m]]$upper.CL <- NaN
  listEffectSizes.MonthTemp.means[[m]]$lower.SEM <- NaN
  listEffectSizes.MonthTemp.means[[m]]$upper.SEM <- NaN
  
  for (idx in 1:nrow(listEffectSizes.MonthTemp.means[[m]])) {
    if (listEffectSizes.MonthTemp.means[[m]]$monthGroup[idx] == refLevel.effectSize.month) {
      # reference
      listEffectSizes.MonthTemp.means[[m]]$mean[idx] <- 1
      listEffectSizes.MonthTemp.means[[m]]$SE[idx] <- listResults.MonthTemp.means[[m]]$SE[idx] / listResults.MonthTemp.means[[m]]$emmean[idx]
      listEffectSizes.MonthTemp.means[[m]]$lower.CL[idx] <- listResults.MonthTemp.means[[m]]$lower.CL[idx] / listResults.MonthTemp.means[[m]]$emmean[idx]
      listEffectSizes.MonthTemp.means[[m]]$upper.CL[idx] <- listResults.MonthTemp.means[[m]]$upper.CL[idx] / listResults.MonthTemp.means[[m]]$emmean[idx]
      listEffectSizes.MonthTemp.means[[m]]$lower.SEM[idx] <- listResults.MonthTemp.means[[m]]$lower.SEM[idx] / listResults.MonthTemp.means[[m]]$emmean[idx]
      listEffectSizes.MonthTemp.means[[m]]$upper.SEM[idx] <- listResults.MonthTemp.means[[m]]$upper.SEM[idx] / listResults.MonthTemp.means[[m]]$emmean[idx]
    } else {
      # other level: calculate uncertainty by linear error term analysis
      # find reference data for same temperature level
      refMask <- (listEffectSizes.MonthTemp.means[[m]]$tempLevel == listEffectSizes.MonthTemp.means[[m]]$tempLevel[idx]) & (listEffectSizes.MonthTemp.means[[m]]$monthGroup == refLevel.effectSize.month)
      # means of reference level and recent data row
      refLevel <- listResults.MonthTemp.means[[m]]$emmean[refMask]
      dataLevel <- listResults.MonthTemp.means[[m]]$emmean[idx]
      
      # error = confidence level
      refDelta <- listResults.MonthTemp.means[[m]]$upper.CL[refMask] - listResults.MonthTemp.means[[m]]$emmean[refMask]
      dataDelta <- listResults.MonthTemp.means[[m]]$upper.CL[idx] - listResults.MonthTemp.means[[m]]$emmean[idx]
      # effect size
      effectSize <- dataLevel / refLevel
      effectSizeDelta <- effectSize * sqrt( (dataDelta / dataLevel)^2 + (refDelta / refLevel)^2 )
      # assign results
      listEffectSizes.MonthTemp.means[[m]]$mean[idx] <- effectSize
      listEffectSizes.MonthTemp.means[[m]]$lower.CL[idx] <- effectSize - effectSizeDelta
      listEffectSizes.MonthTemp.means[[m]]$upper.CL[idx] <- effectSize + effectSizeDelta
      
      # error = SEM
      refDelta <- listResults.MonthTemp.means[[m]]$SE[refMask]
      dataDelta <- listResults.MonthTemp.means[[m]]$SE[idx]
      # effect size
      effectSizeDelta <- effectSize * sqrt( (dataDelta / dataLevel)^2 + (refDelta / refLevel)^2 )
      # assign results
      listEffectSizes.MonthTemp.means[[m]]$SE[idx] <- effectSizeDelta
      listEffectSizes.MonthTemp.means[[m]]$lower.SEM[idx] <- effectSize - effectSizeDelta
      listEffectSizes.MonthTemp.means[[m]]$upper.SEM[idx] <- effectSize + effectSizeDelta
    }
  }
}
rm(refMask, refLevel, dataLevel, refDelta, dataDelta, effectSize, effectSizeDelta)
```


## effect size: fitted values
```{r}
listEffectSizes.fittedValues <- listResults.fittedValues

for (m in 1:length(listModel.final)) {
  listEffectSizes.fittedValues[[m]] <- listEffectSizes.fittedValues[[m]] %>%
    mutate(fitted.normalized = NA)
  
  for (idx in 1:nrow(listEffectSizes.fittedValues[[m]])) {
    # if (listEffectSizes.fittedValues[[m]]$monthGroup[idx] == refLevel.effectSize.month) {
    #   # reference
    #   listEffectSizes.fittedValues[[m]]$fitted.normalized[idx] <- 1
    # } else {
      # other level:
      # find reference data for same animal type
      refMask <- (listResults.MonthTemp.means[[m]]$tempLevel == listEffectSizes.fittedValues[[m]]$tempLevel[idx]) & (listResults.MonthTemp.means[[m]]$monthGroup == refLevel.effectSize.month)
      # means of reference level and recent data row
      refLevel <- listResults.MonthTemp.means[[m]]$emmean[refMask]
      
      listEffectSizes.fittedValues[[m]]$fitted.normalized[idx] <- listEffectSizes.fittedValues[[m]]$fitted[idx] / refLevel
    # }
  }
}

```


## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listEffectSizes.MonthTemp.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listEffectSizes.MonthTemp.contrasts[[m]] <- listEffectSizes.MonthTemp.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value.adj < 0.001 ~ "***",
                               p.value.adj < 0.01  ~ "**",
                               p.value.adj < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value.adj)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1,1,1,2,2,3) + 0.05, times = length(levels.Temp)),
           xEnd     = rep(c(2,3,4,3,4,4) - 0.05, times = length(levels.Temp)),
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots month | temp v2
```{r fig.height=16/2.54, fig.width=10/2.54}
plots.effectSize.MonthTemp.limits <- data.frame(min = c(.8, 0.5, 0.5, 0.8, 0.8, 0.9, 0.5),
                                                max = c(1.3, 2.5, 2.5, 1.75, 1.75, 1.17, 1.7))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listEffectSizes.MonthTemp.contrasts[[m]] <- listEffectSizes.MonthTemp.contrasts[[m]] %>%
    mutate(yStart = rep((0.69 + 0.08*c(0,1,3,0,2,0)) * (plots.effectSize.MonthTemp.limits$max[m] - plots.effectSize.MonthTemp.limits$min[m]) + plots.effectSize.MonthTemp.limits$min[m], times = length(levels.Temp)),
           yEnd = yStart,
           yText = yStart + 0.002*(plots.effectSize.MonthTemp.limits$max[m] - plots.effectSize.MonthTemp.limits$min[m]))
}

plots.effectSize.MonthTemp <- list()
tbl.plotResults.effectSize.MonthTemp <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.effectSize.MonthTemp[[m]] <- ggplot(data = listEffectSizes.MonthTemp.means[[m]])
  
  # define layout parameters
  plots.effectSize.MonthTemp[[m]] <- plots.effectSize.MonthTemp[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 11, angle = 90),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
    facet_grid(rows = vars(tempLevel))
  
  # create point plot: replicate means
  plots.effectSize.MonthTemp[[m]] <- plots.effectSize.MonthTemp[[m]] +
    geom_beeswarm(data = listEffectSizes.fittedValues[[m]],
                  mapping = aes(x = monthGroup,
                                y = fitted.normalized),
                  colour = "#C0C0C0",
                  size = 2,
                  cex = 1.6)
  
  # create LMM plot: mean
  plots.effectSize.MonthTemp[[m]] <- plots.effectSize.MonthTemp[[m]] +
    geom_point(mapping = aes(x = monthGroup,
                             y = mean),
               stat = "identity", # important since we provide mean values already
               size = 3,
               colour = "#1f77bd")
  
  # add LMM plot: error bars
  plots.effectSize.MonthTemp[[m]] <- plots.effectSize.MonthTemp[[m]] +
    geom_errorbar(mapping = aes(x = monthGroup,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = "#1f77bd",
                  position = position_dodge(0.1))
  
  # add discrete x-axis labels
  plots.effectSize.MonthTemp[[m]] <- plots.effectSize.MonthTemp[[m]] +
    scale_x_discrete(labels = labels.Month.Short)

  # add p-values
  plots.effectSize.MonthTemp[[m]] <- plots.effectSize.MonthTemp[[m]] +
    geom_segment(data = listEffectSizes.MonthTemp.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listEffectSizes.MonthTemp.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.effectSize.MonthTemp[[m]] <- plots.effectSize.MonthTemp[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.effectSize.MonthTemp[[m]] <- plots.effectSize.MonthTemp[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.effectSize.MonthTemp.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Month,
         y = listLabelEffectSize[m])
  
  # change coordinates = "zoom in"
  plots.effectSize.MonthTemp[[m]] <- plots.effectSize.MonthTemp[[m]] +
    coord_cartesian(ylim = c(plots.effectSize.MonthTemp.limits$min[m], plots.effectSize.MonthTemp.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_effectSize", "_MonthTemp", ".svg")),
         width = 10, height = 16, units = "cm", dpi = 600, device = "svg",
         plot = plots.effectSize.MonthTemp[[m]], scale = 1)
  
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_effectSize", "_MonthTemp", ".png")),
         width = 10, height = 16, units = "cm", dpi = 600, device = "png",
         plot = plots.effectSize.MonthTemp[[m]], scale = 1)
  
  # results table
  tbl.plotResults.effectSize.MonthTemp <- rbind(tbl.plotResults.effectSize.MonthTemp, data.frame(rep(listOutcome[m], times = nrow(listEffectSizes.MonthTemp.means[[1]])), listEffectSizes.MonthTemp.means[[m]]$monthGroup, listEffectSizes.MonthTemp.means[[m]]$tempLevel, listEffectSizes.MonthTemp.means[[m]]$mean, listEffectSizes.MonthTemp.means[[m]]$SE, listEffectSizes.MonthTemp.means[[m]]$upper.CL - listEffectSizes.MonthTemp.means[[m]]$mean))
}

colnames(tbl.plotResults.effectSize.MonthTemp) <- c("targetVar", "monthGroup", "tempLevel", "mean", "SEM", "CI")
write.csv(tbl.plotResults.effectSize.MonthTemp, file.path(pathWD, paste0("tblPlotResults", "_effectSize", "_MonthTemp", ".csv")), row.names = FALSE)

plots.effectSize.MonthTemp

```


# Characteristic time for month | tempLevel
```{r}
index.tau2 <- which(listOutcome == "tau2")
title.tau2 <- expression("Recovery time" ~ tau["c"] ~ "(ms)")
colorPalette.tau2 <- c("#000000","#1F77BD","#FF7F00","#C21618")
shapePalette.tau2 <- c(21,25,23,24)

listResults.tau2.TempMonth.means <- rbind(listResults.TempMonth.means[[index.tau2]], mutate(listResults.Month.means[[index.tau2]], tempLevel = "Full dataset")) %>%
  mutate(monthGroupNum = case_when(monthGroup == "4-5" ~ 1,
                                   monthGroup == "8"   ~ 2,
                                   monthGroup == "9" ~ 3,
                                   monthGroup == "10" ~ 4,
                                   TRUE ~ NaN)) %>%
  mutate(monthGroupNum = case_when(tempLevel == "Full dataset" ~ monthGroupNum - 0.3,
                                   tempLevel == "Cold" ~ monthGroupNum - 0.1,
                                   tempLevel == "RT" ~ monthGroupNum + 0.1,
                                   tempLevel == "Warm" ~ monthGroupNum + 0.3,
                                   TRUE ~ monthGroupNum))
# sort levels of tempLevel
listResults.tau2.TempMonth.means$tempLevel <- factor(listResults.tau2.TempMonth.means$tempLevel, levels = c("Full dataset","Cold","RT","Warm"))
```

## all plots tau2
```{r fig.height=10/2.54, fig.width=12/2.54}
plots.tau2.TempMonth.limits <- data.frame(min = c(1),
                                          max = c(5))

tbl.plotResults.tau2.TempMonth <- data.frame(matrix(ncol = 5, nrow = 0))

plots.tau2.TempMonth <- ggplot(data = listResults.tau2.TempMonth.means)

# define layout parameters
plots.tau2.TempMonth <- plots.tau2.TempMonth +
  theme_bw() +
  theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
        axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
        axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
        axis.text.y = element_text(size = 10, face = "bold"),
        panel.border = element_rect(size = 1.5, colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
        plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
        plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
        legend.justification = c(1,0),
        legend.position = c(.98,.02),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.text.align = 0)

# create LMM plot: mean
plots.tau2.TempMonth <- plots.tau2.TempMonth +
  geom_point(mapping = aes(x = monthGroupNum,
                           y = emmean,
                           shape = tempLevel),
             stat = "identity", # important since we provide mean values already
             size = 3,
             colour = "#1F77BD",
             fill = "#1F77BD")

# add LMM plot: error bars
plots.tau2.TempMonth <- plots.tau2.TempMonth +
  geom_errorbar(mapping = aes(x = monthGroupNum,
                              ymin = lower.SEM,
                              ymax = upper.SEM),
                size = 0.75,
                width = 0.15,
                colour = "#1F77BD",
                position = position_dodge())

# add discrete x-axis labels
plots.tau2.TempMonth <- plots.tau2.TempMonth +
  scale_x_continuous(breaks = c(1,2,3,4),
                     labels = labels.Month)

# add further layout parameter
plots.tau2.TempMonth <- plots.tau2.TempMonth +
  scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                     limits = c(0, plots.tau2.TempMonth.limits$max),
                     expand = expansion(mult = c(0.025, 0.025))) +
  labs(caption = paste("Model:", format(listFormula.full)),
       x = title.Month,
       y = title.tau2) +
  scale_shape_manual(values = shapePalette.tau2,
                     labels = labels.TempUnit)
  # scale_color_manual(values = colorPalette.tau2,
  #                    labels = labels.TempUnit)

# change coordinates = "zoom in"
plots.tau2.TempMonth <- plots.tau2.TempMonth +
  coord_cartesian(ylim = c(plots.tau2.TempMonth.limits$min, plots.tau2.TempMonth.limits$max))

# save plots
ggsave(filename = file.path(pathWD, paste0(listOutcome[index.tau2], "_v2", "_TempMonth", ".svg")),
       width = 12, height = 10, units = "cm", dpi = 600, device = "svg",
       plot = plots.tau2.TempMonth, scale = 1)

ggsave(filename = file.path(pathWD, paste0(listOutcome[index.tau2], "_v2", "_TempMonth", ".png")),
       width = 12, height = 10, units = "cm", dpi = 600, device = "png",
       plot = plots.tau2.TempMonth, scale = 1)

# results table
tbl.plotResults.tau2.TempMonth <- rbind(tbl.plotResults.tau2.TempMonth, data.frame(rep(listOutcome[index.tau2], times = nrow(listResults.tau2.TempMonth.means)), listResults.tau2.TempMonth.means$tempLevel, listResults.tau2.TempMonth.means$monthGroup, listResults.tau2.TempMonth.means$monthGroupNum, listResults.tau2.TempMonth.means$emmean, listResults.tau2.TempMonth.means$SE))

colnames(tbl.plotResults.tau2.TempMonth) <- c("targetVar", "tempLevel", "monthGroup", "monthGroupNum", "mean", "SEM")
write.csv(tbl.plotResults.tau2.TempMonth, file.path(pathWD, paste0("tblPlotResults", "_v2", "_TempMonth", ".csv")), row.names = FALSE)

plots.tau2.TempMonth
```


# save data space image
```{r}
save.image(file.path(pathWD, paste0("RDataSpace_v9", ".RData")))
```
